diff --git a/applications/solvers/combustion/PDRFoam/Make/options b/applications/solvers/combustion/PDRFoam/Make/options
index 4e5b8fb9a0..a1469d0ef0 100644
--- a/applications/solvers/combustion/PDRFoam/Make/options
+++ b/applications/solvers/combustion/PDRFoam/Make/options
@@ -15,8 +15,7 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -lengine \
diff --git a/applications/test/globalMeshData/Make/files b/applications/test/globalMeshData/Make/files
new file mode 100644
index 0000000000..dc5877e383
--- /dev/null
+++ b/applications/test/globalMeshData/Make/files
@@ -0,0 +1,3 @@
+globalMeshDataTest.C
+
+EXE = $(FOAM_USER_APPBIN)/globalMeshDataTest
diff --git a/applications/test/globalMeshData/Make/options b/applications/test/globalMeshData/Make/options
new file mode 100644
index 0000000000..4c3dd783cb
--- /dev/null
+++ b/applications/test/globalMeshData/Make/options
@@ -0,0 +1,3 @@
+EXE_INC =
+
+EXE_LIBS =
diff --git a/applications/test/globalMeshData/globalMeshDataTest.C b/applications/test/globalMeshData/globalMeshDataTest.C
new file mode 100644
index 0000000000..75f98b1b9c
--- /dev/null
+++ b/applications/test/globalMeshData/globalMeshDataTest.C
@@ -0,0 +1,221 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    globalMeshDataTest
+
+Description
+    Test global point communication
+
+\*---------------------------------------------------------------------------*/
+
+#include "globalMeshData.H"
+#include "argList.H"
+#include "polyMesh.H"
+#include "Time.H"
+#include "mapDistribute.H"
+
+using namespace Foam;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createPolyMesh.H"
+
+    const globalMeshData& globalData = mesh.globalData();
+    const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
+
+
+    // Test:print shared points
+    {
+        const labelListList& globalPointSlaves =
+            globalData.globalPointSlaves();
+        const mapDistribute& globalPointSlavesMap =
+            globalData.globalPointSlavesMap();
+
+        pointField coords(globalPointSlavesMap.constructSize());
+        SubList<point>(coords, coupledPatch.nPoints()).assign
+        (
+            coupledPatch.localPoints()
+        );
+
+        // Exchange data
+        globalPointSlavesMap.distribute(coords);
+
+        // Print
+        forAll(globalPointSlaves, pointI)
+        {
+            const labelList& slavePoints = globalPointSlaves[pointI];
+
+            if (slavePoints.size() > 0)
+            {
+                Pout<< "Master point:" << pointI
+                    << " coord:" << coords[pointI]
+                    << " connected to slave points:" << endl;
+
+                forAll(slavePoints, i)
+                {
+                    Pout<< "    " << coords[slavePoints[i]] << endl;
+                }
+            }
+        }
+    }
+
+
+
+    // Test: point to faces addressing
+    {
+        const labelListList& globalPointBoundaryFaces =
+            globalData.globalPointBoundaryFaces();
+        const mapDistribute& globalPointBoundaryFacesMap =
+            globalData.globalPointBoundaryFacesMap();
+
+        label nBnd = mesh.nFaces()-mesh.nInternalFaces();
+
+        pointField fc(globalPointBoundaryFacesMap.constructSize());
+        SubList<point>(fc, nBnd).assign
+        (
+            primitivePatch
+            (
+                SubList<face>
+                (
+                    mesh.faces(),
+                    nBnd,
+                    mesh.nInternalFaces()
+                ),
+                mesh.points()
+            ).faceCentres()
+        );
+
+        // Exchange data
+        globalPointBoundaryFacesMap.distribute(fc);
+
+        // Print
+        forAll(globalPointBoundaryFaces, pointI)
+        {
+            const labelList& bFaces = globalPointBoundaryFaces[pointI];
+
+            Pout<< "Point:" << pointI
+                << " at:" << coupledPatch.localPoints()[pointI]
+                << " connected to faces:" << endl;
+
+            forAll(bFaces, i)
+            {
+                Pout<< "    " << fc[bFaces[i]] << endl;
+            }
+        }
+    }
+
+
+
+
+
+    // Test:point to cells addressing
+    {
+        const labelList& boundaryCells = globalData.boundaryCells();
+        const labelListList& globalPointBoundaryCells =
+            globalData.globalPointBoundaryCells();
+        const mapDistribute& globalPointBoundaryCellsMap =
+            globalData.globalPointBoundaryCellsMap();
+
+        pointField cc(globalPointBoundaryCellsMap.constructSize());
+        forAll(boundaryCells, i)
+        {
+            cc[i] = mesh.cellCentres()[boundaryCells[i]];
+        }
+
+        // Exchange data
+        globalPointBoundaryCellsMap.distribute(cc);
+
+        // Print
+        forAll(globalPointBoundaryCells, pointI)
+        {
+            const labelList& bCells = globalPointBoundaryCells[pointI];
+
+            Pout<< "Point:" << pointI
+                << " at:" << coupledPatch.localPoints()[pointI]
+                << " connected to cells:" << endl;
+
+            forAll(bCells, i)
+            {
+                Pout<< "    " << cc[bCells[i]] << endl;
+            }
+        }
+    }
+
+
+
+    // Test:print shared edges
+    {
+        const labelListList& globalEdgeSlaves =
+            globalData.globalEdgeSlaves();
+        const mapDistribute& globalEdgeSlavesMap =
+            globalData.globalEdgeSlavesMap();
+
+        // Test: distribute edge centres
+        pointField ec(globalEdgeSlavesMap.constructSize());
+        forAll(coupledPatch.edges(), edgeI)
+        {
+            ec[edgeI] = coupledPatch.edges()[edgeI].centre
+            (
+                coupledPatch.localPoints()
+            );
+        }
+
+        // Exchange data
+        globalEdgeSlavesMap.distribute(ec);
+
+        // Print
+        forAll(globalEdgeSlaves, edgeI)
+        {
+            const labelList& slaveEdges = globalEdgeSlaves[edgeI];
+
+            if (slaveEdges.size() > 0)
+            {
+                Pout<< "Master edge:" << edgeI
+                    << " centre:" << ec[edgeI]
+                    << " connected to slave edges:" << endl;
+
+                forAll(slaveEdges, i)
+                {
+                    Pout<< "    " << ec[slaveEdges[i]] << endl;
+                }
+            }
+        }
+    }
+
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/syncTools/Make/files b/applications/test/syncTools/Make/files
new file mode 100644
index 0000000000..6a186ca84a
--- /dev/null
+++ b/applications/test/syncTools/Make/files
@@ -0,0 +1,3 @@
+syncToolsTest.C
+
+EXE = $(FOAM_USER_APPBIN)/syncToolsTest
diff --git a/applications/test/syncTools/Make/options b/applications/test/syncTools/Make/options
new file mode 100644
index 0000000000..41306609f2
--- /dev/null
+++ b/applications/test/syncTools/Make/options
@@ -0,0 +1 @@
+EXE_INC =
diff --git a/applications/test/syncTools/syncToolsTest.C b/applications/test/syncTools/syncToolsTest.C
new file mode 100644
index 0000000000..320b4e6f61
--- /dev/null
+++ b/applications/test/syncTools/syncToolsTest.C
@@ -0,0 +1,629 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    syncToolsTest
+
+Description
+    Test some functionality in syncTools.
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "argList.H"
+#include "polyMesh.H"
+#include "Time.H"
+#include "Random.H"
+#include "PackedList.H"
+#include "syncTools.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void testPackedList(const polyMesh& mesh, Random& rndGen)
+{
+    Info<< nl << "Testing PackedList synchronisation." << endl;
+
+    {
+        PackedList<3> bits(mesh.nEdges());
+        forAll(bits, i)
+        {
+            bits.set(i, rndGen.integer(0,3));
+        }
+
+        labelList edgeValues(mesh.nEdges());
+        forAll(bits, i)
+        {
+            edgeValues[i] = bits.get(i);
+        }
+
+        PackedList<3> maxBits(bits);
+        labelList maxEdgeValues(edgeValues);
+
+        syncTools::syncEdgeList(mesh, bits, minEqOp<unsigned int>(), 0);
+        syncTools::syncEdgeList(mesh, edgeValues, minEqOp<label>(), 0, false);
+
+        syncTools::syncEdgeList(mesh, maxBits, maxEqOp<unsigned int>(), 0);
+        syncTools::syncEdgeList(mesh, maxEdgeValues, maxEqOp<label>(), 0, false);
+
+        forAll(bits, i)
+        {
+            if
+            (
+                edgeValues[i] != label(bits.get(i))
+             || maxEdgeValues[i] != label(maxBits.get(i))
+            )
+            {
+                FatalErrorIn("testPackedList()")
+                    << "edge:" << i
+                    << " minlabel:" << edgeValues[i]
+                    << " minbits:" << bits.get(i)
+                    << " maxLabel:" << maxEdgeValues[i]
+                    << " maxBits:" << maxBits.get(i)
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    {
+        PackedList<3> bits(mesh.nPoints());
+        forAll(bits, i)
+        {
+            bits.set(i, rndGen.integer(0,3));
+        }
+
+        labelList pointValues(mesh.nPoints());
+        forAll(bits, i)
+        {
+            pointValues[i] = bits.get(i);
+        }
+
+        PackedList<3> maxBits(bits);
+        labelList maxPointValues(pointValues);
+
+        syncTools::syncPointList(mesh, bits, minEqOp<unsigned int>(), 0);
+        syncTools::syncPointList(mesh, pointValues, minEqOp<label>(), 0, false);
+
+        syncTools::syncPointList(mesh, maxBits, maxEqOp<unsigned int>(), 0);
+        syncTools::syncPointList(mesh, maxPointValues, maxEqOp<label>(), 0, false);
+
+        forAll(bits, i)
+        {
+            if
+            (
+                pointValues[i] != label(bits.get(i))
+             || maxPointValues[i] != label(maxBits.get(i))
+            )
+            {
+                FatalErrorIn("testPackedList()")
+                    << "point:" << i
+                    << " minlabel:" << pointValues[i]
+                    << " minbits:" << bits.get(i)
+                    << " maxLabel:" << maxPointValues[i]
+                    << " maxBits:" << maxBits.get(i)
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    {
+        PackedList<3> bits(mesh.nFaces());
+        forAll(bits, faceI)
+        {
+            bits.set(faceI, rndGen.integer(0,3));
+        }
+
+        labelList faceValues(mesh.nFaces());
+        forAll(bits, faceI)
+        {
+            faceValues[faceI] = bits.get(faceI);
+        }
+
+        PackedList<3> maxBits(bits);
+        labelList maxFaceValues(faceValues);
+
+        syncTools::syncFaceList(mesh, bits, minEqOp<unsigned int>());
+        syncTools::syncFaceList(mesh, faceValues, minEqOp<label>(), false);
+
+        syncTools::syncFaceList(mesh, maxBits, maxEqOp<unsigned int>());
+        syncTools::syncFaceList(mesh, maxFaceValues, maxEqOp<label>(), false);
+
+        forAll(bits, faceI)
+        {
+            if
+            (
+                faceValues[faceI] != label(bits.get(faceI))
+             || maxFaceValues[faceI] != label(maxBits.get(faceI))
+            )
+            {
+                FatalErrorIn("testPackedList()")
+                    << "face:" << faceI
+                    << " minlabel:" << faceValues[faceI]
+                    << " minbits:" << bits.get(faceI)
+                    << " maxLabel:" << maxFaceValues[faceI]
+                    << " maxBits:" << maxBits.get(faceI)
+                    << exit(FatalError);
+            }
+        }
+    }
+}
+
+
+void testSparseData(const polyMesh& mesh, Random& rndGen)
+{
+    Info<< nl << "Testing Map synchronisation." << endl;
+
+    primitivePatch allBoundary
+    (
+        SubList<face>
+        (
+            mesh.faces(),
+            mesh.nFaces()-mesh.nInternalFaces(),
+            mesh.nInternalFaces()
+        ),
+        mesh.points()
+    );
+    const pointField& localPoints = allBoundary.localPoints();
+
+    const point greatPoint(GREAT, GREAT, GREAT);
+
+
+    // Point data
+    // ~~~~~~~~~~
+
+    {
+        // Create some data. Use slightly perturbed positions.
+        Map<vector> sparseData;
+        pointField fullData(mesh.nPoints(), greatPoint);
+
+        forAll(localPoints, i)
+        {
+            const point pt = localPoints[i] + 1E-4*rndGen.vector01();
+
+            label meshPointI = allBoundary.meshPoints()[i];
+
+            sparseData.insert(meshPointI, pt);
+            fullData[meshPointI] = pt;
+        }    
+
+        //Pout<< "sparseData:" << sparseData << endl;
+
+        syncTools::syncPointMap
+        (
+            mesh,
+            sparseData,
+            minEqOp<vector>(),
+            true                    // apply separation
+        );
+        syncTools::syncPointList
+        (
+            mesh,
+            fullData,
+            minEqOp<vector>(),
+            greatPoint,
+            true                    // apply separation
+        );
+
+        // Compare.
+        // 1. Is all fullData also present in sparseData and same value
+        forAll(fullData, meshPointI)
+        {
+            const point& fullPt = fullData[meshPointI];
+
+            if (fullPt != greatPoint)
+            {
+                const point& sparsePt = sparseData[meshPointI];
+
+                if (fullPt != sparsePt)
+                {
+                    FatalErrorIn("testSparseData()")
+                        << "point:" << meshPointI
+                        << " full:" << fullPt
+                        << " sparse:" << sparsePt
+                        << exit(FatalError);
+                }
+            }
+        }
+
+        // 2. Does sparseData contain more?
+        forAllConstIter(Map<vector>, sparseData, iter)
+        {
+            const point& sparsePt = iter();
+            label meshPointI = iter.key();
+            const point& fullPt = fullData[meshPointI];
+
+            if (fullPt != sparsePt)
+            {
+                FatalErrorIn("testSparseData()")
+                    << "point:" << meshPointI
+                    << " full:" << fullPt
+                    << " sparse:" << sparsePt
+                    << exit(FatalError);
+            }
+        }
+    }
+
+
+    // Edge data
+    // ~~~~~~~~~
+
+    {
+        // Create some data. Use slightly perturbed positions.
+        EdgeMap<vector> sparseData;
+        pointField fullData(mesh.nEdges(), greatPoint);
+
+        const edgeList& edges = allBoundary.edges();
+        const labelList meshEdges = allBoundary.meshEdges
+        (
+            mesh.edges(),
+            mesh.pointEdges()
+        );
+
+        forAll(edges, i)
+        {
+            const edge& e = edges[i];
+
+            const point pt = e.centre(localPoints) + 1E-4*rndGen.vector01();
+
+            label meshEdgeI = meshEdges[i];
+
+            sparseData.insert(mesh.edges()[meshEdgeI], pt);
+            fullData[meshEdgeI] = pt;
+        }    
+
+        //Pout<< "sparseData:" << sparseData << endl;
+
+        syncTools::syncEdgeMap
+        (
+            mesh,
+            sparseData,
+            minEqOp<vector>(),
+            true
+        );
+        syncTools::syncEdgeList
+        (
+            mesh,
+            fullData,
+            minEqOp<vector>(),
+            greatPoint,
+            true
+        );
+
+        // Compare.
+        // 1. Is all fullData also present in sparseData and same value
+        forAll(fullData, meshEdgeI)
+        {
+            const point& fullPt = fullData[meshEdgeI];
+
+            if (fullPt != greatPoint)
+            {
+                const point& sparsePt = sparseData[mesh.edges()[meshEdgeI]];
+
+                if (fullPt != sparsePt)
+                {
+                    FatalErrorIn("testSparseData()")
+                        << "edge:" << meshEdgeI
+                        << " points:" << mesh.edges()[meshEdgeI]
+                        << " full:" << fullPt
+                        << " sparse:" << sparsePt
+                        << exit(FatalError);
+                }
+            }
+        }
+
+        // 2. Does sparseData contain more?
+        forAll(fullData, meshEdgeI)
+        {
+            const edge& e = mesh.edges()[meshEdgeI];
+
+            EdgeMap<vector>::const_iterator iter = sparseData.find(e);
+
+            if (iter != sparseData.end())
+            {
+                const point& sparsePt = iter();
+                const point& fullPt = fullData[meshEdgeI];
+
+                if (fullPt != sparsePt)
+                {
+                    FatalErrorIn("testSparseData()")
+                        << "Extra edge:" << meshEdgeI
+                        << " points:" << mesh.edges()[meshEdgeI]
+                        << " full:" << fullPt
+                        << " sparse:" << sparsePt
+                        << exit(FatalError);
+                }
+            }
+        }
+    }
+}
+
+
+void testPointSync(const polyMesh& mesh, Random& rndGen)
+{
+    Info<< nl << "Testing point-wise data synchronisation." << endl;
+
+    const point greatPoint(GREAT, GREAT, GREAT);
+
+    // Test position.
+
+    {
+        WarningIn("testPointSync()")
+            << "Position test only correct for cases without cyclics"
+            << " with shared points." << endl;
+
+        pointField syncedPoints(mesh.points());
+        syncTools::syncPointList
+        (
+            mesh,
+            syncedPoints,
+            minEqOp<point>(),
+            greatPoint,
+            true
+        );
+
+        forAll(syncedPoints, pointI)
+        {
+            if (mag(syncedPoints[pointI] - mesh.points()[pointI]) > SMALL)
+            {
+                FatalErrorIn("testPointSync()")
+                    << "Point " << pointI
+                    << " original location " << mesh.points()[pointI]
+                    << " synced location " << syncedPoints[pointI]
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    // Test masterPoints
+
+    {
+        labelList nMasters(mesh.nPoints(), 0);
+
+        PackedBoolList isMasterPoint(syncTools::getMasterPoints(mesh));
+
+        forAll(isMasterPoint, pointI)
+        {
+            if (isMasterPoint[pointI])
+            {
+                nMasters[pointI] = 1;
+            }
+        }
+
+        syncTools::syncPointList
+        (
+            mesh,
+            nMasters,
+            plusEqOp<label>(),
+            0,
+            false
+        );
+
+        forAll(nMasters, pointI)
+        {
+            if (nMasters[pointI] != 1)
+            {
+                //FatalErrorIn("testPointSync()")
+                WarningIn("testPointSync()")
+                    << "Point " << pointI
+                    << " original location " << mesh.points()[pointI]
+                    << " has " << nMasters[pointI]
+                    << " masters."
+                    //<< exit(FatalError);
+                    << endl;
+            }
+        }
+    }
+}
+
+
+void testEdgeSync(const polyMesh& mesh, Random& rndGen)
+{
+    Info<< nl << "Testing edge-wise data synchronisation." << endl;
+
+    const edgeList& edges = mesh.edges();
+
+    const point greatPoint(GREAT, GREAT, GREAT);
+
+    // Test position.
+
+    {
+        WarningIn("testEdgeSync()")
+            << "Position test only correct for cases without cyclics"
+            << " with shared edges." << endl;
+
+        pointField syncedMids(edges.size());
+        forAll(syncedMids, edgeI)
+        {
+            syncedMids[edgeI] = edges[edgeI].centre(mesh.points());
+        }
+        syncTools::syncEdgeList
+        (
+            mesh,
+            syncedMids,
+            minEqOp<point>(),
+            greatPoint,
+            true
+        );
+
+        forAll(syncedMids, edgeI)
+        {
+            point eMid = edges[edgeI].centre(mesh.points());
+
+            if (mag(syncedMids[edgeI] - eMid) > SMALL)
+            {
+                FatalErrorIn("testEdgeSync()")
+                    << "Edge " << edgeI
+                    << " original midpoint " << eMid
+                    << " synced location " << syncedMids[edgeI]
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    // Test masterEdges
+
+    {
+        labelList nMasters(edges.size(), 0);
+
+        PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+
+        forAll(isMasterEdge, edgeI)
+        {
+            if (isMasterEdge[edgeI])
+            {
+                nMasters[edgeI] = 1;
+            }
+        }
+
+        syncTools::syncEdgeList
+        (
+            mesh,
+            nMasters,
+            plusEqOp<label>(),
+            0,
+            false
+        );
+
+        forAll(nMasters, edgeI)
+        {
+            if (nMasters[edgeI] != 1)
+            {
+                //FatalErrorIn("testEdgeSync()")
+                WarningIn("testEdgeSync()")
+                    << "Edge " << edgeI
+                    << " midpoint " << edges[edgeI].centre(mesh.points())
+                    << " has " << nMasters[edgeI]
+                    << " masters."
+                    //<< exit(FatalError);
+                    << endl;
+            }
+        }
+    }
+}
+
+
+void testFaceSync(const polyMesh& mesh, Random& rndGen)
+{
+    Info<< nl << "Testing face-wise data synchronisation." << endl;
+
+    // Test position.
+
+    {
+        pointField syncedFc(mesh.faceCentres());
+
+        syncTools::syncFaceList
+        (
+            mesh,
+            syncedFc,
+            maxEqOp<point>(),
+            true
+        );
+
+        forAll(syncedFc, faceI)
+        {
+            if (mag(syncedFc[faceI] - mesh.faceCentres()[faceI]) > SMALL)
+            {
+                FatalErrorIn("testFaceSync()")
+                    << "Face " << faceI
+                    << " original centre " << mesh.faceCentres()[faceI]
+                    << " synced centre " << syncedFc[faceI]
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    // Test masterFaces
+
+    {
+        labelList nMasters(mesh.nFaces(), 0);
+
+        PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh));
+
+        forAll(isMasterFace, faceI)
+        {
+            if (isMasterFace[faceI])
+            {
+                nMasters[faceI] = 1;
+            }
+        }
+
+        syncTools::syncFaceList
+        (
+            mesh,
+            nMasters,
+            plusEqOp<label>(),
+            false
+        );
+
+        forAll(nMasters, faceI)
+        {
+            if (nMasters[faceI] != 1)
+            {
+                FatalErrorIn("testFaceSync()")
+                    << "Face " << faceI
+                    << " centre " << mesh.faceCentres()[faceI]
+                    << " has " << nMasters[faceI]
+                    << " masters."
+                    << exit(FatalError);
+            }
+        }
+    }
+}
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createPolyMesh.H"
+
+
+    Random rndGen(5341*(Pstream::myProcNo()+1));
+
+
+    // Face sync
+    testFaceSync(mesh, rndGen);
+
+    // Edge sync
+    testEdgeSync(mesh, rndGen);
+
+    // Point sync
+    testPointSync(mesh, rndGen);
+
+    // PackedList synchronisation
+    testPackedList(mesh, rndGen);
+
+    // Sparse synchronisation
+    testSparseData(mesh, rndGen);
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/advanced/modifyMesh/Make/options b/applications/utilities/mesh/advanced/modifyMesh/Make/options
index 39929ff890..70c838b774 100644
--- a/applications/utilities/mesh/advanced/modifyMesh/Make/options
+++ b/applications/utilities/mesh/advanced/modifyMesh/Make/options
@@ -1,9 +1,7 @@
 EXE_INC = \
-    -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
-    -ltriSurface \
     -lmeshTools \
     -ldynamicMesh
diff --git a/applications/utilities/mesh/advanced/splitCells/Make/options b/applications/utilities/mesh/advanced/splitCells/Make/options
index ca424359cf..7349856cab 100644
--- a/applications/utilities/mesh/advanced/splitCells/Make/options
+++ b/applications/utilities/mesh/advanced/splitCells/Make/options
@@ -1,9 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
     -ldynamicMesh \
-    -ltriSurface \
     -lmeshTools
diff --git a/applications/utilities/mesh/conversion/foamToSurface/Make/files b/applications/utilities/mesh/conversion/foamToSurface/Make/files
new file mode 100644
index 0000000000..a5dfd1e846
--- /dev/null
+++ b/applications/utilities/mesh/conversion/foamToSurface/Make/files
@@ -0,0 +1,3 @@
+foamToSurface.C
+
+EXE = $(FOAM_APPBIN)/foamToSurface
diff --git a/applications/utilities/mesh/conversion/foamToSurface/Make/options b/applications/utilities/mesh/conversion/foamToSurface/Make/options
new file mode 100644
index 0000000000..a504dd8617
--- /dev/null
+++ b/applications/utilities/mesh/conversion/foamToSurface/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/surfMesh/lnInclude
+
+EXE_LIBS = \
+    -lsurfMesh
diff --git a/applications/utilities/mesh/conversion/foamToSurface/foamToSurface.C b/applications/utilities/mesh/conversion/foamToSurface/foamToSurface.C
new file mode 100644
index 0000000000..f9a255c502
--- /dev/null
+++ b/applications/utilities/mesh/conversion/foamToSurface/foamToSurface.C
@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2010-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    foamToSurface
+
+Description
+    Reads an OpenFOAM mesh and writes the boundaries in a surface format.
+
+Usage
+    - foamToSurface [OPTION] \n
+    Reads an OpenFOAM mesh and writes the boundaries in a surface format.
+
+    @param -scale \<factor\>\n
+    Specify an alternative geometry scaling factor.
+    Eg, use @b 1000 to scale @em [m] to @em [mm].
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "timeSelector.H"
+#include "Time.H"
+#include "polyMesh.H"
+
+#include "MeshedSurfaces.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noParallel();
+    argList::validArgs.append("outputFile.ext");
+    timeSelector::addOptions();
+
+    argList::addOption
+    (
+        "scale",
+        "scale",
+        "specify geometry scaling factor"
+    );
+
+#   include "setRootCase.H"
+
+    const stringList& params = args.additionalArgs();
+
+    scalar scaleFactor = 0;
+    args.optionReadIfPresent<scalar>("scale", scaleFactor);
+
+    fileName exportName(params[0]);
+
+    fileName exportBase = exportName.lessExt();
+    word exportExt = exportName.ext();
+
+    if (!meshedSurface::canWriteType(exportExt, true))
+    {
+        return 1;
+    }
+
+#   include "createTime.H"
+    instantList timeDirs = timeSelector::select0(runTime, args);
+#   include "createPolyMesh.H"
+
+    forAll(timeDirs, timeI)
+    {
+        runTime.setTime(timeDirs[timeI], timeI);
+#       include "getTimeIndex.H"
+
+        polyMesh::readUpdateState state = mesh.readUpdate();
+
+        if (timeI == 0 || state != polyMesh::UNCHANGED)
+        {
+            if (state == polyMesh::UNCHANGED)
+            {
+                exportName = exportBase + "." + exportExt;
+            }
+            else
+            {
+                exportName =
+                    exportBase + '_' + runTime.timeName() + "." + exportExt;
+            }
+
+            meshedSurface surf(mesh.boundaryMesh());
+            surf.scalePoints(scaleFactor);
+
+            Info<< "writing " << exportName;
+            if (scaleFactor <= 0)
+            {
+                Info<< " without scaling" << endl;
+            }
+            else
+            {
+                Info<< " with scaling " << scaleFactor << endl;
+            }
+            surf.write(exportName);
+        }
+
+        Info<< nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/conversion/foamToSurface/getTimeIndex.H b/applications/utilities/mesh/conversion/foamToSurface/getTimeIndex.H
new file mode 100644
index 0000000000..85a92405dd
--- /dev/null
+++ b/applications/utilities/mesh/conversion/foamToSurface/getTimeIndex.H
@@ -0,0 +1,51 @@
+// Read time index from */uniform/time, but treat 0 and constant specially
+
+    word timeName = "0";
+
+    if
+    (
+        runTime.timeName() != "constant"
+     && runTime.timeName() != "0"
+    )
+    {
+        IOobject io
+        (
+            "time",
+            runTime.timeName(),
+            "uniform",
+            runTime,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        if (io.headerOk())
+        {
+            IOdictionary timeObject
+            (
+                IOobject
+                (
+                    "time",
+                    runTime.timeName(),
+                    "uniform",
+                    runTime,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE,
+                    false
+                )
+            );
+
+            label index;
+            timeObject.lookup("index") >> index;
+            timeName = Foam::name(index);
+        }
+        else
+        {
+            timeName = runTime.timeName();
+            // Info<< "skip ... missing entry " << io.objectPath() << endl;
+            // continue;
+        }
+    }
+
+    Info<< "\nTime [" << timeName << "] = " << runTime.timeName() << nl;
+
diff --git a/applications/utilities/mesh/conversion/ideasUnvToFoam/Make/options b/applications/utilities/mesh/conversion/ideasUnvToFoam/Make/options
index 9f08e8d2a8..02c293ceed 100644
--- a/applications/utilities/mesh/conversion/ideasUnvToFoam/Make/options
+++ b/applications/utilities/mesh/conversion/ideasUnvToFoam/Make/options
@@ -1,7 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \
-    -ltriSurface
+    -lsurfMesh
diff --git a/applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C b/applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C
index ba7ee39c80..259ccd2069 100644
--- a/applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C
+++ b/applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C
@@ -41,7 +41,8 @@ Description
 #include "cellSet.H"
 #include "faceSet.H"
 #include "DynamicList.H"
-#include "triSurface.H"
+
+#include "MeshedSurfaces.H"
 
 using namespace Foam;
 
@@ -268,9 +269,11 @@ void readCells
             break;
         }
 
-        IStringStream lineStr(line);
         label cellI, feID, physProp, matProp, colour, nNodes;
-        lineStr >> cellI >> feID >> physProp >> matProp >> colour >> nNodes;
+
+        IStringStream lineStr(line);
+        lineStr
+            >> cellI >> feID >> physProp >> matProp >> colour >> nNodes;
 
         if (foundFeType.insert(feID))
         {
@@ -297,7 +300,8 @@ void readCells
 
             face cVerts(3);
             IStringStream lineStr(line);
-            lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2];
+            lineStr
+                >> cVerts[0] >> cVerts[1] >> cVerts[2];
             boundaryFaces.append(cVerts);
             boundaryFaceIndices.append(cellI);
         }
@@ -308,7 +312,8 @@ void readCells
 
             face cVerts(4);
             IStringStream lineStr(line);
-            lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
+            lineStr
+                >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
             boundaryFaces.append(cVerts);
             boundaryFaceIndices.append(cellI);
         }
@@ -319,14 +324,15 @@ void readCells
 
             labelList cVerts(4);
             IStringStream lineStr(line);
-            lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
+            lineStr
+                >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3];
 
             cellVerts.append(cellShape(tet, cVerts, true));
             cellMaterial.append(physProp);
 
             if (cellVerts.last().size() != cVerts.size())
             {
-                Pout<< "Line:" << is.lineNumber()
+                Info<< "Line:" << is.lineNumber()
                     << " element:" << cellI
                     << " type:" << feID
                     << " collapsed from " << cVerts << nl
@@ -341,15 +347,16 @@ void readCells
 
             labelList cVerts(6);
             IStringStream lineStr(line);
-            lineStr >> cVerts[0] >> cVerts[1] >> cVerts[2] >> cVerts[3]
-                    >> cVerts[4] >> cVerts[5];
+            lineStr
+                >> cVerts[0] >> cVerts[1] >> cVerts[2]
+                >> cVerts[3] >> cVerts[4] >> cVerts[5];
 
             cellVerts.append(cellShape(prism, cVerts, true));
             cellMaterial.append(physProp);
 
             if (cellVerts.last().size() != cVerts.size())
             {
-                Pout<< "Line:" << is.lineNumber()
+                Info<< "Line:" << is.lineNumber()
                     << " element:" << cellI
                     << " type:" << feID
                     << " collapsed from " << cVerts << nl
@@ -373,7 +380,7 @@ void readCells
 
             if (cellVerts.last().size() != cVerts.size())
             {
-                Pout<< "Line:" << is.lineNumber()
+                Info<< "Line:" << is.lineNumber()
                     << " element:" << cellI
                     << " type:" << feID
                     << " collapsed from " << cVerts << nl
@@ -388,7 +395,7 @@ void readCells
                 IOWarningIn("readCells(IFstream&, label&)", is)
                     << "Cell type " << feID << " not supported" << endl;
             }
-            is.getLine(line);  //Do nothing
+            is.getLine(line);  // Do nothing
         }
     }
 
@@ -579,7 +586,11 @@ int main(int argc, char *argv[])
 {
     argList::noParallel();
     argList::validArgs.append(".unv file");
-    argList::addBoolOption("dump");
+    argList::addBoolOption
+    (
+        "dump",
+        "dump boundary faces as boundaryFaces.obj (for debugging)"
+    );
 
 #   include "setRootCase.H"
 #   include "createTime.H"
@@ -858,40 +869,25 @@ int main(int argc, char *argv[])
     polyPoints /= lengthScale;
 
 
-    // For debugging: dump boundary faces as triSurface
+    // For debugging: dump boundary faces as OBJ surface mesh
     if (args.optionFound("dump"))
     {
-        DynamicList<labelledTri> triangles(boundaryFaces.size());
-
-        forAll(boundaryFaces, i)
-        {
-            const face& f = boundaryFaces[i];
-
-            faceList triFaces(f.nTriangles(polyPoints));
-            label nTri = 0;
-            f.triangles(polyPoints, nTri, triFaces);
-
-            forAll(triFaces, triFaceI)
-            {
-                const face& f = triFaces[triFaceI];
-                triangles.append(labelledTri(f[0], f[1], f[2], 0));
-            }
-        }
-
-        // Create globally numbered tri surface
-        triSurface rawSurface(triangles.shrink(), polyPoints);
-
-        // Create locally numbered tri surface
-        triSurface surface
-        (
-            rawSurface.localFaces(),
-            rawSurface.localPoints()
-        );
-
-        Info<< "Writing boundary faces to STL file boundaryFaces.stl"
+        Info<< "Writing boundary faces to OBJ file boundaryFaces.obj"
             << nl << endl;
 
-        surface.write(runTime.path()/"boundaryFaces.stl");
+        // Create globally numbered surface
+        meshedSurface rawSurface
+        (
+            xferCopy(polyPoints),
+            xferCopyTo< faceList >(boundaryFaces)
+        );
+
+        // Write locally numbered surface
+        meshedSurface
+        (
+            xferCopy(rawSurface.localPoints()),
+            xferCopy(rawSurface.localFaces())
+        ).write(runTime.path()/"boundaryFaces.obj");
     }
 
 
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index ab8adca478..c7f2d5b6f3 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -180,6 +180,12 @@ castellatedMeshControls
     // NOTE: This point should never be on a face, always inside a cell, even
     // after refinement.
     locationInMesh (5 0.28 0.43);
+
+
+    // Whether any faceZones (as specified in the refinementSurfaces)
+    // are only on the boundary of corresponding cellZones or also allow
+    // free-standing zone faces. Not used if there are no faceZones.
+    allowFreeStandingZoneFaces true;
 }
 
 
@@ -274,6 +280,7 @@ addLayersControls
     maxThicknessToMedialRatio 0.3;
 
     // Angle used to pick up medial axis points
+    // Note: changed(corrected) w.r.t 16x! 90 degrees corresponds to 130 in 16x.
     minMedianAxisAngle 90;
 
     // Create buffer region for new layer terminations
diff --git a/applications/utilities/mesh/manipulation/autoPatch/Make/options b/applications/utilities/mesh/manipulation/autoPatch/Make/options
index d002df8030..7349856cab 100644
--- a/applications/utilities/mesh/manipulation/autoPatch/Make/options
+++ b/applications/utilities/mesh/manipulation/autoPatch/Make/options
@@ -1,10 +1,7 @@
 EXE_INC = \
-    -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
-
 EXE_LIBS = \
-    -ltriSurface \
     -ldynamicMesh \
     -lmeshTools
diff --git a/applications/utilities/mesh/manipulation/objToVTK/Make/options b/applications/utilities/mesh/manipulation/objToVTK/Make/options
index 0548bdd022..79a9626b55 100644
--- a/applications/utilities/mesh/manipulation/objToVTK/Make/options
+++ b/applications/utilities/mesh/manipulation/objToVTK/Make/options
@@ -1,8 +1 @@
-/*
-EXE_INC = \
-    -I$(LIB_SRC)/triSurface/lnInclude
-
-EXE_LIBS = \
-    -lmeshTools \
-    -ltriSurface
-*/
+/* */
diff --git a/applications/utilities/mesh/manipulation/refineMesh/Make/options b/applications/utilities/mesh/manipulation/refineMesh/Make/options
index 5c9c1c2b87..73f34b0f85 100644
--- a/applications/utilities/mesh/manipulation/refineMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/refineMesh/Make/options
@@ -1,9 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
-
 EXE_LIBS = \
     -ldynamicMesh \
     -lmeshTools
diff --git a/applications/utilities/miscellaneous/foamDebugSwitches/Make/options b/applications/utilities/miscellaneous/foamDebugSwitches/Make/options
index e473b84bb9..7207b53927 100644
--- a/applications/utilities/miscellaneous/foamDebugSwitches/Make/options
+++ b/applications/utilities/miscellaneous/foamDebugSwitches/Make/options
@@ -18,6 +18,7 @@ EXE_LIBS = \
     -lengine \
     -lerrorEstimation \
     -lfieldFunctionObjects \
+    -lfileFormats \
     -lfiniteVolume \
     -lforces \
     -lfvMotionSolvers \
diff --git a/applications/utilities/postProcessing/sampling/probeLocations/Make/options b/applications/utilities/postProcessing/sampling/probeLocations/Make/options
index 44392ffc44..ec8bc31179 100644
--- a/applications/utilities/postProcessing/sampling/probeLocations/Make/options
+++ b/applications/utilities/postProcessing/sampling/probeLocations/Make/options
@@ -2,7 +2,6 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 EXE_LIBS = \
@@ -10,5 +9,4 @@ EXE_LIBS = \
     -lgenericPatchFields \
     -lmeshTools \
     -lsampling \
-    -ltriSurface \
     -llagrangian
diff --git a/applications/utilities/postProcessing/sampling/sample/Make/options b/applications/utilities/postProcessing/sampling/sample/Make/options
index 4877a53003..88c6039b4e 100644
--- a/applications/utilities/postProcessing/sampling/sample/Make/options
+++ b/applications/utilities/postProcessing/sampling/sample/Make/options
@@ -3,7 +3,6 @@ EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude
 
 EXE_LIBS = \
@@ -12,5 +11,4 @@ EXE_LIBS = \
     -lmeshTools \
     -lsampling \
     -lsurfMesh \
-    -ltriSurface \
     -llagrangian
diff --git a/applications/utilities/surface/surfaceAutoPatch/Make/options b/applications/utilities/surface/surfaceAutoPatch/Make/options
index e4751082e1..9f08e8d2a8 100644
--- a/applications/utilities/surface/surfaceAutoPatch/Make/options
+++ b/applications/utilities/surface/surfaceAutoPatch/Make/options
@@ -1,5 +1,4 @@
 EXE_INC = \
-    /* -I$(LIB_SRC)/cfdTools/general/lnInclude */ \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude
 
diff --git a/applications/utilities/surface/surfaceFeatureConvert/Make/options b/applications/utilities/surface/surfaceFeatureConvert/Make/options
index f87dafaa5d..e71834f99e 100644
--- a/applications/utilities/surface/surfaceFeatureConvert/Make/options
+++ b/applications/utilities/surface/surfaceFeatureConvert/Make/options
@@ -1,9 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/edgeMesh/lnInclude \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/edgeMesh/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \
-    -ledgeMesh \
-    -ltriSurface
+    -ledgeMesh
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index f38bf13cc5..71eda554fa 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -79,11 +79,7 @@ void deleteBox
     {
         const point eMid = surf.edges()[edgeI].centre(surf.localPoints());
 
-        if
-        (
-            (removeInside && bb.contains(eMid))
-         || (!removeInside && !bb.contains(eMid))
-        )
+        if (removeInside ? bb.contains(eMid) : !bb.contains(eMid))
         {
             edgeStat[edgeI] = surfaceFeatures::NONE;
         }
@@ -133,7 +129,6 @@ int main(int argc, char *argv[])
 
 
 
-
     // Either construct features from surface&featureangle or read set.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
diff --git a/applications/utilities/surface/surfaceFind/Make/options b/applications/utilities/surface/surfaceFind/Make/options
index 4c41bc1a49..946b16a9ad 100644
--- a/applications/utilities/surface/surfaceFind/Make/options
+++ b/applications/utilities/surface/surfaceFind/Make/options
@@ -1,6 +1,6 @@
 EXE_INC = \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \
-    -ltriSurface
+    -lsurfMesh
diff --git a/applications/utilities/surface/surfaceFind/surfaceFind.C b/applications/utilities/surface/surfaceFind/surfaceFind.C
index dcaccc04c2..8e36fc25ec 100644
--- a/applications/utilities/surface/surfaceFind/surfaceFind.C
+++ b/applications/utilities/surface/surfaceFind/surfaceFind.C
@@ -23,18 +23,16 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Finds nearest triangle and vertex.
+    Finds nearest face and vertex.
 
 \*---------------------------------------------------------------------------*/
 
-#include "triSurface.H"
 #include "argList.H"
 #include "OFstream.H"
 
-#ifndef namespaceFoam
-#define namespaceFoam
-    using namespace Foam;
-#endif
+#include "MeshedSurfaces.H"
+
+using namespace Foam;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,9 +42,9 @@ int main(int argc, char *argv[])
 {
     argList::noParallel();
     argList::validArgs.clear();
-    argList::addOption("x", "X");
-    argList::addOption("y", "Y");
-    argList::addOption("z", "Z");
+    argList::addOption("x", "X", "The point x-coordinate (if non-zero)");
+    argList::addOption("y", "Y", "The point y-coordinate (if non-zero)");
+    argList::addOption("z", "Z", "The point y-coordinate (if non-zero)");
 
     argList::validArgs.append("surface file");
 
@@ -54,15 +52,15 @@ int main(int argc, char *argv[])
 
     point samplePt
     (
-        args.optionRead<scalar>("x"),
-        args.optionRead<scalar>("y"),
-        args.optionRead<scalar>("z")
+        args.optionLookupOrDefault<scalar>("x", 0),
+        args.optionLookupOrDefault<scalar>("y", 0),
+        args.optionLookupOrDefault<scalar>("z", 0)
     );
     Info<< "Looking for nearest face/vertex to " << samplePt << endl;
 
 
-    Info<< "Reading surf1 ..." << endl;
-    triSurface surf1(args.additionalArgs()[0]);
+    Info<< "Reading surf ..." << endl;
+    meshedSurface surf1(args.additionalArgs()[0]);
 
     //
     // Nearest vertex
@@ -83,11 +81,11 @@ int main(int argc, char *argv[])
         }
     }
 
-    Info<< "Nearest vertex:" << endl
-        << "    index      :" << minIndex << " (in localPoints)" << endl
+    Info<< "Nearest vertex:" << nl
+        << "    index      :" << minIndex << " (in localPoints)" << nl
         << "    index      :" << surf1.meshPoints()[minIndex]
-        << " (in points)" << endl
-        << "    coordinates:" << localPoints[minIndex] << endl
+        << " (in points)" << nl
+        << "    coordinates:" << localPoints[minIndex] << nl
         << endl;
 
     //
@@ -101,8 +99,7 @@ int main(int argc, char *argv[])
 
     forAll(surf1, faceI)
     {
-        const labelledTri& f = surf1[faceI];
-        const point centre = f.centre(points);
+        const point centre = surf1[faceI].centre(points);
 
         const scalar dist = mag(centre - samplePt);
         if (dist < minDist)
@@ -112,16 +109,19 @@ int main(int argc, char *argv[])
         }
     }
 
-    const labelledTri& f = surf1[minIndex];
+    const face& f = surf1[minIndex];
 
-    Info<< "Face with nearest centre:" << endl
-        << "    index        :" << minIndex << endl
-        << "    centre       :" << f.centre(points) << endl
-        << "    face         :" << f << endl
-        << "    vertex coords:" << points[f[0]] << " "
-        << points[f[1]] << " " << points[f[2]] << endl
-        << endl;
+    Info<< "Face with nearest centre:" << nl
+        << "    index        :" << minIndex << nl
+        << "    centre       :" << f.centre(points) << nl
+        << "    face         :" << f << nl
+        << "    vertex coords:\n";
+    forAll(f, fp)
+    {
+        Info<< "        " << points[f[fp]] << "\n";
+    }
 
+    Info<< endl;
 
     Info<< "End\n" << endl;
 
diff --git a/applications/utilities/surface/surfaceInertia/surfaceInertia.C b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
index 12785fec9b..f3bcb4599b 100644
--- a/applications/utilities/surface/surfaceInertia/surfaceInertia.C
+++ b/applications/utilities/surface/surfaceInertia/surfaceInertia.C
@@ -202,6 +202,7 @@ void massPropertiesSolid
     J *= density;
 }
 
+
 void massPropertiesShell
 (
     const pointField& pts,
@@ -272,7 +273,7 @@ tensor applyParallelAxisTheorem
 
     vector d = (refPt - cM);
 
-    return (J + m*((d & d)*I - d*d));
+    return J + m*((d & d)*I - d*d);
 }
 
 
@@ -389,7 +390,7 @@ int main(int argc, char *argv[])
         << eVec.x() <<  ' ' << eVec.y() << ' ' << eVec.z()
         << endl;
 
-    if(calcAroundRefPt)
+    if (calcAroundRefPt)
     {
         Info << "Inertia tensor relative to " << refPt << " = "
             << applyParallelAxisTheorem(m, cM, J, refPt)
diff --git a/applications/utilities/surface/surfaceSmooth/Make/options b/applications/utilities/surface/surfaceSmooth/Make/options
index 4c41bc1a49..946b16a9ad 100644
--- a/applications/utilities/surface/surfaceSmooth/Make/options
+++ b/applications/utilities/surface/surfaceSmooth/Make/options
@@ -1,6 +1,6 @@
 EXE_INC = \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 EXE_LIBS = \
     -lmeshTools \
-    -ltriSurface
+    -lsurfMesh
diff --git a/applications/utilities/surface/surfaceSmooth/surfaceSmooth.C b/applications/utilities/surface/surfaceSmooth/surfaceSmooth.C
index bce23c9846..e68f6ab20b 100644
--- a/applications/utilities/surface/surfaceSmooth/surfaceSmooth.C
+++ b/applications/utilities/surface/surfaceSmooth/surfaceSmooth.C
@@ -23,15 +23,16 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Example of simple laplacian smoother
+    Example of a simple laplacian smoother
 
 \*---------------------------------------------------------------------------*/
 
-#include "triSurface.H"
 #include "argList.H"
 #include "OFstream.H"
 #include "boundBox.H"
 
+#include "MeshedSurfaces.H"
+
 using namespace Foam;
 
 
@@ -61,23 +62,20 @@ int main(int argc, char *argv[])
     label iters(readLabel(IStringStream(args.additionalArgs()[2])()));
     fileName outFileName(args.additionalArgs()[3]);
 
-    Info<< "Relax:" << relax << endl;
-    Info<< "Iters:" << iters << endl;
+    Info<< "Relax:" << relax << nl
+        << "Iters:" << iters << nl
+        << "Reading surface from " << surfFileName << " ..." << endl;
 
+    meshedSurface surf1(surfFileName);
 
-    Info<< "Reading surface from " << surfFileName << " ..." << endl;
-
-    triSurface surf1(surfFileName);
-
-    Info<< "Triangles    : " << surf1.size() << endl;
-    Info<< "Vertices     : " << surf1.nPoints() << endl;
-    Info<< "Bounding Box : " << boundBox(surf1.localPoints()) << endl;
+    Info<< "Faces        : " << surf1.size() << nl
+        << "Vertices     : " << surf1.nPoints() << nl
+        << "Bounding Box : " << boundBox(surf1.localPoints()) << endl;
 
     pointField newPoints(surf1.localPoints());
 
     const labelListList& pointEdges = surf1.pointEdges();
 
-
     for (label iter = 0; iter < iters; iter++)
     {
         forAll(pointEdges, vertI)
@@ -100,16 +98,14 @@ int main(int argc, char *argv[])
         }
     }
 
-    triSurface surf2
-    (
-        surf1.localFaces(),
-        surf1.patches(),
-        newPoints
-    );
-
     Info<< "Writing surface to " << outFileName << " ..." << endl;
 
-    surf2.write(outFileName);
+    meshedSurface
+    (
+        xferMove(newPoints),
+        xferCopy(surf1.localFaces()),
+        xferCopy(surf1.surfZones())
+    ).write(outFileName);
 
     Info<< "End\n" << endl;
 
diff --git a/applications/utilities/surface/surfaceSubset/surfaceSubset.C b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
index dd8a6ede00..6cc8ee4678 100644
--- a/applications/utilities/surface/surfaceSubset/surfaceSubset.C
+++ b/applications/utilities/surface/surfaceSubset/surfaceSubset.C
@@ -100,10 +100,8 @@ int main(int argc, char *argv[])
         meshSubsetDict.lookup("addFaceNeighbours")
     );
 
-    Switch invertSelection
-    (
-        meshSubsetDict.lookup("invertSelection")
-    );
+    const bool invertSelection =
+        meshSubsetDict.lookupOrDefault<bool>("invertSelection", false);
 
     // Mark the cells for the subset
 
@@ -246,7 +244,7 @@ int main(int argc, char *argv[])
         // bb of surface
         treeBoundBox bb(selectSurf.localPoints());
 
-        // Radnom number generator
+        // Random number generator
         Random rndGen(354543);
 
         // search engine
@@ -269,14 +267,11 @@ int main(int argc, char *argv[])
                 indexedOctree<treeDataTriSurface>::volumeType t =
                     selectTree.getVolumeType(fc);
 
-                if (t == indexedOctree<treeDataTriSurface>::INSIDE && !outside)
-                {
-                    facesToSubset[faceI] = true;
-                }
-                else if
+                if
                 (
-                    t == indexedOctree<treeDataTriSurface>::OUTSIDE
-                 && outside
+                    outside
+                  ? (t == indexedOctree<treeDataTriSurface>::OUTSIDE)
+                  : (t == indexedOctree<treeDataTriSurface>::INSIDE)
                 )
                 {
                     facesToSubset[faceI] = true;
@@ -346,20 +341,11 @@ int main(int argc, char *argv[])
     if (invertSelection)
     {
         Info<< "Inverting selection." << endl;
-        boolList newFacesToSubset(facesToSubset.size());
 
         forAll(facesToSubset, i)
         {
-            if (facesToSubset[i])
-            {
-                newFacesToSubset[i] = false;
-            }
-            else
-            {
-                newFacesToSubset[i] = true;
-            }
+            facesToSubset[i] = facesToSubset[i] ? false : true;
         }
-        facesToSubset.transfer(newFacesToSubset);
     }
 
 
diff --git a/applications/utilities/surface/surfaceSubset/surfaceSubsetDict b/applications/utilities/surface/surfaceSubset/surfaceSubsetDict
index 3a8db2cffd..f956a24789 100644
--- a/applications/utilities/surface/surfaceSubset/surfaceSubsetDict
+++ b/applications/utilities/surface/surfaceSubset/surfaceSubsetDict
@@ -36,7 +36,7 @@ zone
 surface
 {
     name "sphere.stl";
-    outside yes;
+    outside     yes;
 }
 
 // Extend selection with edge neighbours
diff --git a/applications/utilities/surface/surfaceTransformPoints/Make/options b/applications/utilities/surface/surfaceTransformPoints/Make/options
index 4c41bc1a49..a504dd8617 100644
--- a/applications/utilities/surface/surfaceTransformPoints/Make/options
+++ b/applications/utilities/surface/surfaceTransformPoints/Make/options
@@ -1,6 +1,5 @@
 EXE_INC = \
-    -I$(LIB_SRC)/triSurface/lnInclude
+    -I$(LIB_SRC)/surfMesh/lnInclude
 
 EXE_LIBS = \
-    -lmeshTools \
-    -ltriSurface
+    -lsurfMesh
diff --git a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
index 8bb87d3ee9..387471521e 100644
--- a/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
+++ b/applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C
@@ -35,7 +35,6 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "triSurface.H"
 #include "argList.H"
 #include "OFstream.H"
 #include "IFstream.H"
@@ -45,6 +44,8 @@ Description
 #include "quaternion.H"
 #include "mathematicalConstants.H"
 
+#include "MeshedSurfaces.H"
+
 using namespace Foam;
 using namespace Foam::constant::mathematical;
 
@@ -114,7 +115,7 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
-    triSurface surf1(surfFileName);
+    meshedSurface surf1(surfFileName);
 
     pointField points(surf1.points());
 
@@ -193,9 +194,8 @@ int main(int argc, char *argv[])
         points.replace(vector::Z, scaleVector.z()*points.component(vector::Z));
     }
 
-    triSurface surf2(surf1, surf1.patches(), points);
-
-    surf2.write(outFileName);
+    surf1.movePoints(points);
+    surf1.write(outFileName);
 
     Info<< "End\n" << endl;
 
diff --git a/src/Allwmake b/src/Allwmake
index 969687eb83..f24f2f370d 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -22,12 +22,14 @@ wmake libso OpenFOAM
 
 wmake libso lagrangian/basic
 
+wmake libso fileFormats
 wmake libso edgeMesh
 wmake libso surfMesh
 wmake libso triSurface
 
 # Decomposition methods needed by meshTools
 wmake libso parallel/decompositionMethods
+wmake libso parallel/metisDecomp
 
 wmake libso meshTools
 wmake libso finiteVolume
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index 147b370d1a..1b3bb264a2 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -32,6 +32,8 @@ License
 #include "demandDrivenData.H"
 #include "globalPoints.H"
 //#include "geomGlobalPoints.H"
+#include "polyMesh.H"
+#include "mapDistribute.H"
 #include "labelIOList.H"
 #include "PackedList.H"
 #include "mergePoints.H"
@@ -115,24 +117,17 @@ void Foam::globalMeshData::initProcAddr()
 // Given information about locally used edges allocate global shared edges.
 void Foam::globalMeshData::countSharedEdges
 (
-    const HashTable<labelList, edge, Hash<edge> >& procSharedEdges,
-    HashTable<label, edge, Hash<edge> >& globalShared,
+    const EdgeMap<labelList>& procSharedEdges,
+    EdgeMap<label>& globalShared,
     label& sharedEdgeI
 )
 {
     // Count occurrences of procSharedEdges in global shared edges table.
-    for
-    (
-        HashTable<labelList, edge, Hash<edge> >::const_iterator iter =
-            procSharedEdges.begin();
-        iter != procSharedEdges.end();
-        ++iter
-    )
+    forAllConstIter(EdgeMap<labelList>, procSharedEdges, iter)
     {
         const edge& e = iter.key();
 
-        HashTable<label, edge, Hash<edge> >::iterator globalFnd =
-            globalShared.find(e);
+        EdgeMap<label>::iterator globalFnd = globalShared.find(e);
 
         if (globalFnd == globalShared.end())
         {
@@ -189,10 +184,7 @@ void Foam::globalMeshData::calcSharedEdges() const
     // Find edges using shared points. Store correspondence to local edge
     // numbering. Note that multiple local edges can have the same shared
     // points! (for cyclics or separated processor patches)
-    HashTable<labelList, edge, Hash<edge> > localShared
-    (
-        2*sharedPtAddr.size()
-    );
+    EdgeMap<labelList> localShared(2*sharedPtAddr.size());
 
     const edgeList& edges = mesh_.edges();
 
@@ -218,7 +210,7 @@ void Foam::globalMeshData::calcSharedEdges() const
                     sharedPtAddr[e1Fnd()]
                 );
 
-                HashTable<labelList, edge, Hash<edge> >::iterator iter =
+                EdgeMap<labelList>::iterator iter =
                     localShared.find(sharedEdge);
 
                 if (iter == localShared.end())
@@ -249,7 +241,7 @@ void Foam::globalMeshData::calcSharedEdges() const
     // used). But then this only gets done once so not too bothered about the
     // extra global communication.
 
-    HashTable<label, edge, Hash<edge> > globalShared(nGlobalPoints());
+    EdgeMap<label> globalShared(nGlobalPoints());
 
     if (Pstream::master())
     {
@@ -275,10 +267,7 @@ void Foam::globalMeshData::calcSharedEdges() const
             {
                 // Receive the edges using shared points from the slave.
                 IPstream fromSlave(Pstream::blocking, slave);
-                HashTable<labelList, edge, Hash<edge> > procSharedEdges
-                (
-                    fromSlave
-                );
+                EdgeMap<labelList> procSharedEdges(fromSlave);
 
                 if (debug)
                 {
@@ -295,17 +284,11 @@ void Foam::globalMeshData::calcSharedEdges() const
         // These were only used once so are not proper shared edges.
         // Remove them.
         {
-            HashTable<label, edge, Hash<edge> > oldSharedEdges(globalShared);
+            EdgeMap<label> oldSharedEdges(globalShared);
 
             globalShared.clear();
 
-            for
-            (
-                HashTable<label, edge, Hash<edge> >::const_iterator iter =
-                    oldSharedEdges.begin();
-                iter != oldSharedEdges.end();
-                ++iter
-            )
+            forAllConstIter(EdgeMap<label>, oldSharedEdges, iter)
             {
                 if (iter() != -1)
                 {
@@ -361,18 +344,11 @@ void Foam::globalMeshData::calcSharedEdges() const
     DynamicList<label> dynSharedEdgeLabels(globalShared.size());
     DynamicList<label> dynSharedEdgeAddr(globalShared.size());
 
-    for
-    (
-        HashTable<labelList, edge, Hash<edge> >::const_iterator iter =
-            localShared.begin();
-        iter != localShared.end();
-        ++iter
-    )
+    forAllConstIter(EdgeMap<labelList>, localShared, iter)
     {
         const edge& e = iter.key();
 
-        HashTable<label, edge, Hash<edge> >::const_iterator edgeFnd =
-            globalShared.find(e);
+        EdgeMap<label>::const_iterator edgeFnd = globalShared.find(e);
 
         if (edgeFnd != globalShared.end())
         {
@@ -434,6 +410,557 @@ Foam::label Foam::globalMeshData::countCoincidentFaces
 }
 
 
+void Foam::globalMeshData::calcGlobalPointSlaves() const
+{
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalPointSlaves() :"
+            << " calculating coupled master to slave point addressing."
+            << endl;
+    }
+
+    // Calculate connected points for master points
+    globalPoints globalData(mesh_, coupledPatch(), true);
+
+    const Map<label>& meshToProcPoint = globalData.meshToProcPoint();
+
+    // Create global numbering for coupled points
+    globalPointNumberingPtr_.reset
+    (
+        new globalIndex(globalData.globalIndices())
+    );
+    const globalIndex& globalIndices = globalPointNumberingPtr_();
+
+    // Create master to slave addressing. Empty for slave points.
+    globalPointSlavesPtr_.reset
+    (
+        new labelListList(coupledPatch().nPoints())
+    );
+    labelListList& globalPointSlaves = globalPointSlavesPtr_();
+
+    forAllConstIter(Map<label>, meshToProcPoint, iter)
+    {
+        label localPointI = iter.key();
+        const labelList& pPoints = globalData.procPoints()[iter()];
+
+        // Am I master?
+        if
+        (
+            globalIndices.isLocal(pPoints[0])
+         && globalIndices.toLocal(pPoints[0]) == localPointI
+        )
+        {
+            labelList& slaves = globalPointSlaves[localPointI];
+            slaves.setSize(pPoints.size()-1);
+            for (label i = 1; i < pPoints.size(); i++)
+            {
+                slaves[i-1] = pPoints[i];
+            }
+        }
+    }
+
+    // Create schedule to get information from slaves onto master
+
+    // Construct compact numbering and distribution map.
+    // Changes globalPointSlaves to be indices into compact data
+
+    List<Map<label> > compactMap(Pstream::nProcs());
+    globalPointSlavesMapPtr_.reset
+    (
+        new mapDistribute
+        (
+            globalIndices,
+            globalPointSlaves,
+            compactMap
+        )
+    );
+
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalPointSlaves() :"
+            << " coupled points:" << coupledPatch().nPoints()
+            << " additional remote points:"
+            <<  globalPointSlavesMapPtr_().constructSize()
+              - coupledPatch().nPoints()
+            << endl;
+    }
+}
+
+
+void Foam::globalMeshData::calcGlobalEdgeSlaves() const
+{
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
+            << " calculating coupled master to slave edge addressing."
+            << endl;
+    }
+
+    const labelListList& globalPointSlaves = this->globalPointSlaves();
+    const mapDistribute& globalPointSlavesMap = this->globalPointSlavesMap();
+
+    // - Send across connected edges (in global edge addressing)
+    // - Check on receiving side whether edge has same slave edge
+    //   on both endpoints.
+
+    // Create global numbering for coupled edges
+    globalEdgeNumberingPtr_.reset
+    (
+        new globalIndex(coupledPatch().nEdges())
+    );
+    const globalIndex& globalIndices = globalEdgeNumberingPtr_();
+
+    // Coupled point to global coupled edges.
+    labelListList globalPointEdges(globalPointSlavesMap.constructSize());
+
+    // Create local version
+    const labelListList& pointEdges = coupledPatch().pointEdges();
+    forAll(pointEdges, pointI)
+    {
+        const labelList& pEdges = pointEdges[pointI];
+        labelList& globalPEdges = globalPointEdges[pointI];
+        globalPEdges.setSize(pEdges.size());
+        forAll(pEdges, i)
+        {
+            globalPEdges[i] = globalIndices.toGlobal(pEdges[i]);
+        }
+    }
+
+    // Pull slave data to master
+    globalPointSlavesMap.distribute(globalPointEdges);
+
+    // Now check on master if any of my edges are also on slave.
+    // This assumes that if slaves have a coupled edge it is also on
+    // the master (otherwise the mesh would be illegal anyway)
+
+    labelHashSet pointEdgeSet;
+
+    const edgeList& edges = coupledPatch().edges();
+
+    // Create master to slave addressing. Empty for slave edges.
+    globalEdgeSlavesPtr_.reset(new labelListList(edges.size()));
+    labelListList& globalEdgeSlaves = globalEdgeSlavesPtr_();
+
+    forAll(edges, edgeI)
+    {
+        const edge& e = edges[edgeI];
+        const labelList& slaves0 = globalPointSlaves[e[0]];
+        const labelList& slaves1 = globalPointSlaves[e[1]];
+
+        // Check for edges that are in both slaves0 and slaves1.
+        pointEdgeSet.clear();
+        forAll(slaves0, i)
+        {
+            const labelList& connectedEdges = globalPointEdges[slaves0[i]];
+            pointEdgeSet.insert(connectedEdges);
+        }
+        forAll(slaves1, i)
+        {
+            const labelList& connectedEdges = globalPointEdges[slaves1[i]];
+
+            forAll(connectedEdges, j)
+            {
+                label globalEdgeI = connectedEdges[j];
+
+                if (pointEdgeSet.found(globalEdgeI))
+                {
+                    // Found slave edge.
+                    label sz = globalEdgeSlaves[edgeI].size();
+                    globalEdgeSlaves[edgeI].setSize(sz+1);
+                    globalEdgeSlaves[edgeI][sz] = globalEdgeI;
+                }
+            }
+        }
+    }
+
+
+    // Construct map
+    List<Map<label> > compactMap(Pstream::nProcs());
+    globalEdgeSlavesMapPtr_.reset
+    (
+        new mapDistribute
+        (
+            globalIndices,
+            globalEdgeSlaves,
+            compactMap
+        )
+    );
+
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalEdgeSlaves() :"
+            << " coupled edge:" << edges.size()
+            << " additional remote edges:"
+            << globalEdgeSlavesMapPtr_().constructSize() - edges.size()
+            << endl;
+    }
+}
+
+
+// Calculate uncoupled boundary faces (without calculating
+// primitiveMesh::pointFaces())
+void Foam::globalMeshData::calcPointBoundaryFaces
+(
+    labelListList& pointBoundaryFaces
+) const
+{
+    const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
+    const Map<label>& meshPointMap = coupledPatch().meshPointMap();
+
+    // 1. Count
+
+    labelList nPointFaces(coupledPatch().nPoints(), 0);
+
+    forAll(bMesh, patchI)
+    {
+        const polyPatch& pp = bMesh[patchI];
+
+        if (!pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                const face& f = pp[i];
+
+                forAll(f, fp)
+                {
+                    Map<label>::const_iterator iter = meshPointMap.find(f[fp]);
+                    if (iter != meshPointMap.end())
+                    {
+                        nPointFaces[iter()]++;
+                    }
+                }
+            }
+        }
+    }
+
+
+    // 2. Size
+
+    pointBoundaryFaces.setSize(coupledPatch().nPoints());
+    forAll(nPointFaces, pointI)
+    {
+        pointBoundaryFaces[pointI].setSize(nPointFaces[pointI]);
+    }
+    nPointFaces = 0;
+
+
+    // 3. Fill
+
+    forAll(bMesh, patchI)
+    {
+        const polyPatch& pp = bMesh[patchI];
+
+        if (!pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                const face& f = pp[i];
+                forAll(f, fp)
+                {
+                    Map<label>::const_iterator iter = meshPointMap.find(f[fp]);
+                    if (iter != meshPointMap.end())
+                    {
+                        label bFaceI = pp.start() + i - mesh_.nInternalFaces();
+                        pointBoundaryFaces[iter()][nPointFaces[iter()]++] =
+                            bFaceI;
+                    }
+                }
+            }
+        }
+    }
+}
+
+void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
+{
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalPointBoundaryFaces() :"
+            << " calculating coupled point to boundary face addressing."
+            << endl;
+    }
+
+    // Construct local point to (uncoupled)boundaryfaces.
+    labelListList pointBoundaryFaces;
+    calcPointBoundaryFaces(pointBoundaryFaces);
+
+
+    // Global indices for boundary faces
+    globalBoundaryFaceNumberingPtr_.reset
+    (
+        new globalIndex(mesh_.nFaces()-mesh_.nInternalFaces())
+    );
+    globalIndex& globalIndices = globalBoundaryFaceNumberingPtr_();
+
+
+    // Convert local boundary faces to global numbering
+    globalPointBoundaryFacesPtr_.reset
+    (
+        new labelListList(globalPointSlavesMap().constructSize())
+    );
+    labelListList& globalPointBoundaryFaces = globalPointBoundaryFacesPtr_();
+
+    forAll(pointBoundaryFaces, pointI)
+    {
+        const labelList& bFaces = pointBoundaryFaces[pointI];
+        labelList& globalFaces = globalPointBoundaryFaces[pointI];
+        globalFaces.setSize(bFaces.size());
+        forAll(bFaces, i)
+        {
+            globalFaces[i] = globalIndices.toGlobal(bFaces[i]);
+        }
+    }
+
+
+    // Pull slave pointBoundaryFaces to master
+    globalPointSlavesMap().distribute(globalPointBoundaryFaces);
+
+
+    // Merge slave labels into master globalPointBoundaryFaces
+    const labelListList& pointSlaves = globalPointSlaves();
+
+    forAll(pointSlaves, pointI)
+    {
+        const labelList& slaves = pointSlaves[pointI];
+
+        if (slaves.size() > 0)
+        {
+            labelList& myBFaces = globalPointBoundaryFaces[pointI];
+
+            forAll(slaves, i)
+            {
+                const labelList& slaveBFaces =
+                    globalPointBoundaryFaces[slaves[i]];
+
+                // Add all slaveBFaces. Note that need to check for
+                // uniqueness only in case of cyclics.
+
+                label sz = myBFaces.size();
+                myBFaces.setSize(sz+slaveBFaces.size());
+                forAll(slaveBFaces, j)
+                {
+                    label slave = slaveBFaces[j];
+                    if (findIndex(SubList<label>(myBFaces, sz), slave) == -1)
+                    {
+                        myBFaces[sz++] = slave;
+                    }
+                }
+                myBFaces.setSize(sz);
+            }
+        }
+    }
+
+
+    // Copy merged boundaryFaces back from master into slave slot
+    forAll(pointSlaves, pointI)
+    {
+        const labelList& bFaces = globalPointBoundaryFaces[pointI];
+        const labelList& slaves = pointSlaves[pointI];
+
+        forAll(slaves, i)
+        {
+            globalPointBoundaryFaces[slaves[i]] = bFaces;
+        }
+    }
+
+
+    // Sync back to slaves.
+    globalPointSlavesMap().reverseDistribute
+    (
+        coupledPatch().nPoints(),
+        globalPointBoundaryFaces
+    );
+
+
+    // Construct a map to get the face data directly
+    List<Map<label> > compactMap(Pstream::nProcs());
+    globalPointBoundaryFacesMapPtr_.reset
+    (
+        new mapDistribute
+        (
+            globalIndices,
+            globalPointBoundaryFaces,
+            compactMap
+        )
+    );
+
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalPointBoundaryFaces() :"
+            << " coupled points:" << coupledPatch().nPoints()
+            << " local boundary faces:" <<  globalIndices.localSize()
+            << " additional remote faces:"
+            <<  globalPointBoundaryFacesMapPtr_().constructSize()
+              - globalIndices.localSize()
+            << endl;
+    }
+}
+
+
+void Foam::globalMeshData::calcGlobalPointBoundaryCells() const
+{
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalPointBoundaryCells() :"
+            << " calculating coupled point to boundary cell addressing."
+            << endl;
+    }
+
+    // Create map of boundary cells and point-cell addressing
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    label bCellI = 0;
+    Map<label> meshCellMap(4*coupledPatch().nPoints());
+    DynamicList<label> cellMap(meshCellMap.size());
+
+    // Create addressing for point to boundary cells (local)
+    labelListList pointBoundaryCells(coupledPatch().nPoints());
+
+    forAll(coupledPatch().meshPoints(), pointI)
+    {
+        label meshPointI = coupledPatch().meshPoints()[pointI];
+        const labelList& pCells = mesh_.pointCells(meshPointI);
+
+        labelList& bCells = pointBoundaryCells[pointI];
+        bCells.setSize(pCells.size());
+
+        forAll(pCells, i)
+        {
+            label cellI = pCells[i];
+            Map<label>::iterator fnd = meshCellMap.find(cellI);
+
+            if (fnd != meshCellMap.end())
+            {
+                bCells[i] = fnd();
+            }
+            else
+            {
+                meshCellMap.insert(cellI, bCellI);
+                cellMap.append(cellI);
+                bCells[i] = bCellI;
+                bCellI++;
+            }
+        }
+    }
+
+
+    boundaryCellsPtr_.reset(new labelList());
+    labelList& boundaryCells = boundaryCellsPtr_();
+    boundaryCells.transfer(cellMap.shrink());
+
+
+    // Convert point-cells to global point numbers
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    globalBoundaryCellNumberingPtr_.reset
+    (
+        new globalIndex(boundaryCells.size())
+    );
+    globalIndex& globalIndices = globalBoundaryCellNumberingPtr_();
+
+
+    globalPointBoundaryCellsPtr_.reset
+    (
+        new labelListList(globalPointSlavesMap().constructSize())
+    );
+    labelListList& globalPointBoundaryCells = globalPointBoundaryCellsPtr_();
+
+    forAll(pointBoundaryCells, pointI)
+    {
+        const labelList& pCells = pointBoundaryCells[pointI];
+        labelList& globalCells = globalPointBoundaryCells[pointI];
+        globalCells.setSize(pCells.size());
+        forAll(pCells, i)
+        {
+            globalCells[i] = globalIndices.toGlobal(pCells[i]);
+        }
+    }
+
+
+    // Pull slave pointBoundaryCells to master
+    globalPointSlavesMap().distribute(globalPointBoundaryCells);
+
+
+    // Merge slave labels into master globalPointBoundaryCells
+    const labelListList& pointSlaves = globalPointSlaves();
+
+    forAll(pointSlaves, pointI)
+    {
+        const labelList& slaves = pointSlaves[pointI];
+
+        if (slaves.size() > 0)
+        {
+            labelList& myBCells = globalPointBoundaryCells[pointI];
+
+            forAll(slaves, i)
+            {
+                const labelList& slaveBCells =
+                    globalPointBoundaryCells[slaves[i]];
+
+                // Add all slaveBCells. Note that need to check for
+                // uniqueness only in case of cyclics.
+
+                label sz = myBCells.size();
+                myBCells.setSize(sz+slaveBCells.size());
+                forAll(slaveBCells, j)
+                {
+                    label slave = slaveBCells[j];
+                    if (findIndex(SubList<label>(myBCells, sz), slave) == -1)
+                    {
+                        myBCells[sz++] = slave;
+                    }
+                }
+                myBCells.setSize(sz);
+            }
+        }
+    }
+
+
+    // Copy merged boundaryCells back from master into slave slot
+    forAll(pointSlaves, pointI)
+    {
+        const labelList& bCells = globalPointBoundaryCells[pointI];
+        const labelList& slaves = pointSlaves[pointI];
+
+        forAll(slaves, i)
+        {
+            globalPointBoundaryCells[slaves[i]] = bCells;
+        }
+    }
+
+
+    // Sync back to slaves.
+    globalPointSlavesMap().reverseDistribute
+    (
+        coupledPatch().nPoints(),
+        globalPointBoundaryCells
+    );
+
+
+    // Construct a map to get the cell data directly
+    List<Map<label> > compactMap(Pstream::nProcs());
+    globalPointBoundaryCellsMapPtr_.reset
+    (
+        new mapDistribute
+        (
+            globalIndices,
+            globalPointBoundaryCells,
+            compactMap
+        )
+    );
+
+    if (debug)
+    {
+        Pout<< "globalMeshData::calcGlobalPointBoundaryCells() :"
+            << " coupled points:" << coupledPatch().nPoints()
+            << " local boundary cells:" <<  globalIndices.localSize()
+            << " additional remote cells:"
+            <<  globalPointBoundaryCellsMapPtr_().constructSize()
+              - globalIndices.localSize()
+            << endl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from polyMesh
@@ -511,6 +1038,25 @@ void Foam::globalMeshData::clearOut()
     nGlobalPoints_ = -1;
     deleteDemandDrivenData(sharedEdgeLabelsPtr_);
     deleteDemandDrivenData(sharedEdgeAddrPtr_);
+
+    coupledPatchPtr_.clear();
+    // Point
+    globalPointNumberingPtr_.clear();
+    globalPointSlavesPtr_.clear();
+    globalPointSlavesMapPtr_.clear();
+    // Edge
+    globalEdgeNumberingPtr_.clear();
+    globalEdgeSlavesPtr_.clear();
+    globalEdgeSlavesMapPtr_.clear();
+    // Face
+    globalBoundaryFaceNumberingPtr_.clear();
+    globalPointBoundaryFacesPtr_.clear();
+    globalPointBoundaryFacesMapPtr_.clear();
+    // Cell
+    boundaryCellsPtr_.clear();
+    globalBoundaryCellNumberingPtr_.clear();
+    globalPointBoundaryCellsPtr_.clear();
+    globalPointBoundaryCellsMapPtr_.clear();
 }
 
 
@@ -709,6 +1255,203 @@ const Foam::labelList& Foam::globalMeshData::sharedEdgeAddr() const
 }
 
 
+const Foam::indirectPrimitivePatch& Foam::globalMeshData::coupledPatch() const
+{
+    if (!coupledPatchPtr_.valid())
+    {
+        const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
+
+        label nCoupled = 0;
+
+        forAll(bMesh, patchI)
+        {
+            const polyPatch& pp = bMesh[patchI];
+
+            if (pp.coupled())
+            {
+                nCoupled += pp.size();
+            }
+        }
+        labelList coupledFaces(nCoupled);
+        nCoupled = 0;
+
+        forAll(bMesh, patchI)
+        {
+            const polyPatch& pp = bMesh[patchI];
+
+            if (pp.coupled())
+            {
+                label faceI = pp.start();
+
+                forAll(pp, i)
+                {
+                    coupledFaces[nCoupled++] = faceI++;
+                }
+            }
+        }
+
+        coupledPatchPtr_.reset
+        (
+            new indirectPrimitivePatch
+            (
+                IndirectList<face>
+                (
+                    mesh_.faces(),
+                    coupledFaces
+                ),
+                mesh_.points()
+            )
+        );
+
+        if (debug)
+        {
+            Pout<< "globalMeshData::coupledPatch() :"
+                << " constructed  coupled faces patch:"
+                << "  faces:" << coupledPatchPtr_().size()
+                << "  points:" << coupledPatchPtr_().nPoints()
+                << endl;
+        }
+    }
+    return coupledPatchPtr_();
+}
+
+
+const Foam::globalIndex& Foam::globalMeshData::globalPointNumbering() const
+{
+    if (!globalPointNumberingPtr_.valid())
+    {
+        calcGlobalPointSlaves();
+    }
+    return globalPointNumberingPtr_();
+}
+
+
+const Foam::labelListList& Foam::globalMeshData::globalPointSlaves() const
+{
+    if (!globalPointSlavesPtr_.valid())
+    {
+        calcGlobalPointSlaves();
+    }
+    return globalPointSlavesPtr_();
+}
+
+
+const Foam::mapDistribute& Foam::globalMeshData::globalPointSlavesMap() const
+{
+    if (!globalPointSlavesMapPtr_.valid())
+    {
+        calcGlobalPointSlaves();
+    }
+    return globalPointSlavesMapPtr_();
+}
+
+
+const Foam::globalIndex& Foam::globalMeshData::globalEdgeNumbering() const
+{
+    if (!globalEdgeNumberingPtr_.valid())
+    {
+        calcGlobalEdgeSlaves();
+    }
+    return globalEdgeNumberingPtr_();
+}
+
+
+const Foam::labelListList& Foam::globalMeshData::globalEdgeSlaves() const
+{
+    if (!globalEdgeSlavesPtr_.valid())
+    {
+        calcGlobalEdgeSlaves();
+    }
+    return globalEdgeSlavesPtr_();
+}
+
+
+const Foam::mapDistribute& Foam::globalMeshData::globalEdgeSlavesMap() const
+{
+    if (!globalEdgeSlavesMapPtr_.valid())
+    {
+        calcGlobalEdgeSlaves();
+    }
+    return globalEdgeSlavesMapPtr_();
+}
+
+
+const Foam::globalIndex& Foam::globalMeshData::globalBoundaryFaceNumbering()
+const
+{
+    if (!globalBoundaryFaceNumberingPtr_.valid())
+    {
+        calcGlobalPointBoundaryFaces();
+    }
+    return globalBoundaryFaceNumberingPtr_();
+}
+
+
+const Foam::labelListList& Foam::globalMeshData::globalPointBoundaryFaces()
+const
+{
+    if (!globalPointBoundaryFacesPtr_.valid())
+    {
+        calcGlobalPointBoundaryFaces();
+    }
+    return globalPointBoundaryFacesPtr_();
+}
+
+
+const Foam::mapDistribute& Foam::globalMeshData::globalPointBoundaryFacesMap()
+const
+{
+    if (!globalPointBoundaryFacesMapPtr_.valid())
+    {
+        calcGlobalPointBoundaryFaces();
+    }
+    return globalPointBoundaryFacesMapPtr_();
+}
+
+
+const Foam::labelList& Foam::globalMeshData::boundaryCells() const
+{
+    if (!boundaryCellsPtr_.valid())
+    {
+        calcGlobalPointBoundaryCells();
+    }
+    return boundaryCellsPtr_();
+}
+
+
+const Foam::globalIndex& Foam::globalMeshData::globalBoundaryCellNumbering()
+const
+{
+    if (!globalBoundaryCellNumberingPtr_.valid())
+    {
+        calcGlobalPointBoundaryCells();
+    }
+    return globalBoundaryCellNumberingPtr_();
+}
+
+
+const Foam::labelListList& Foam::globalMeshData::globalPointBoundaryCells()
+const
+{
+    if (!globalPointBoundaryCellsPtr_.valid())
+    {
+        calcGlobalPointBoundaryCells();
+    }
+    return globalPointBoundaryCellsPtr_();
+}
+
+
+const Foam::mapDistribute& Foam::globalMeshData::globalPointBoundaryCellsMap()
+const
+{
+    if (!globalPointBoundaryCellsMapPtr_.valid())
+    {
+        calcGlobalPointBoundaryCells();
+    }
+    return globalPointBoundaryCellsMapPtr_();
+}
+
+
 void Foam::globalMeshData::movePoints(const pointField& newPoints)
 {
     // Topology does not change and we don't store any geometry so nothing
@@ -740,7 +1483,7 @@ void Foam::globalMeshData::updateMesh()
     // Option 1. Topological
     {
         // Calculate all shared points. This does all the hard work.
-        globalPoints parallelPoints(mesh_);
+        globalPoints parallelPoints(mesh_, false);
 
         // Copy data out.
         nGlobalPoints_ = parallelPoints.nGlobalPoints();
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
index 7c4893126d..a6e03baf95 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
@@ -74,10 +74,13 @@ SourceFiles
 #ifndef globalMeshData_H
 #define globalMeshData_H
 
-#include "polyMesh.H"
+//#include "polyMesh.H"
 #include "Switch.H"
 #include "processorTopology.H"
 #include "labelPair.H"
+#include "indirectPrimitivePatch.H"
+#include "boundBox.H"
+#include "IOobject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -88,7 +91,10 @@ namespace Foam
 
 class globalMeshData;
 Ostream& operator<<(Ostream&, const globalMeshData&);
-
+class globalIndex;
+class polyMesh;
+class mapDistribute;
+template<class T> class EdgeMap;
 
 /*---------------------------------------------------------------------------*\
                            Class globalMeshData Declaration
@@ -194,6 +200,39 @@ class globalMeshData
             mutable labelList* sharedEdgeAddrPtr_;
 
 
+        // Coupled point addressing
+        // This is addressing from coupled point to coupled points,faces,cells
+        // This is a full schedule so includes points only used by two
+        // coupled patches.
+
+            mutable autoPtr<indirectPrimitivePatch> coupledPatchPtr_;
+
+            // Coupled point to coupled points
+
+                mutable autoPtr<globalIndex> globalPointNumberingPtr_;
+                mutable autoPtr<labelListList> globalPointSlavesPtr_;
+                mutable autoPtr<mapDistribute> globalPointSlavesMapPtr_;
+
+            // Coupled edge to coupled edges
+
+                mutable autoPtr<globalIndex> globalEdgeNumberingPtr_;
+                mutable autoPtr<labelListList> globalEdgeSlavesPtr_;
+                mutable autoPtr<mapDistribute> globalEdgeSlavesMapPtr_;
+
+            // Coupled point to boundary faces
+
+                mutable autoPtr<globalIndex> globalBoundaryFaceNumberingPtr_;
+                mutable autoPtr<labelListList> globalPointBoundaryFacesPtr_;
+                mutable autoPtr<mapDistribute> globalPointBoundaryFacesMapPtr_;
+
+            // Coupled point to boundary cells
+
+                mutable autoPtr<labelList> boundaryCellsPtr_;
+                mutable autoPtr<globalIndex> globalBoundaryCellNumberingPtr_;
+                mutable autoPtr<labelListList> globalPointBoundaryCellsPtr_;
+                mutable autoPtr<mapDistribute> globalPointBoundaryCellsMapPtr_;
+
+
     // Private Member Functions
 
         //- Set up processor patch addressing
@@ -202,8 +241,8 @@ class globalMeshData
         //- Helper function for shared edge addressing
         static void countSharedEdges
         (
-            const HashTable<labelList, edge, Hash<edge> >& procSharedEdges,
-            HashTable<label, edge, Hash<edge> >&,
+            const EdgeMap<labelList>&,
+            EdgeMap<label>&,
             label&
         );
 
@@ -217,6 +256,22 @@ class globalMeshData
             const vectorField& separationDist
         );
 
+        //- Calculate global point addressing.
+        void calcGlobalPointSlaves() const;
+
+        //- Calculate global edge addressing.
+        void calcGlobalEdgeSlaves() const;
+
+        //- Calculate coupled point to uncoupled boundary faces. Local only.
+        void calcPointBoundaryFaces(labelListList& pointBoundaryFaces) const;
+
+        //- Calculate global point to global boundary face addressing.
+        void calcGlobalPointBoundaryFaces() const;
+
+        //- Calculate global point to global boundary cell addressing.
+        void calcGlobalPointBoundaryCells() const;
+
+
         //- Disallow default bitwise copy construct
         globalMeshData(const globalMeshData&);
 
@@ -388,6 +443,47 @@ public:
             const labelList& sharedEdgeAddr() const;
 
 
+
+        // Global master - slave point communication
+
+            //- Return patch of all coupled faces
+            const indirectPrimitivePatch& coupledPatch() const;
+
+            // Coupled point to coupled points. Coupled points are points on
+            // any coupled patch.
+
+                //- Numbering of coupled points is according to coupledPatch.
+                const globalIndex& globalPointNumbering() const;
+                //- For every coupled point the indices into the field
+                //  distributed by below map.
+                const labelListList& globalPointSlaves() const;
+                const mapDistribute& globalPointSlavesMap() const;
+
+            // Coupled edge to coupled edges.
+
+                const globalIndex& globalEdgeNumbering() const;
+                const labelListList& globalEdgeSlaves() const;
+                const mapDistribute& globalEdgeSlavesMap() const;
+
+            // Coupled point to boundary faces. These are uncoupled boundary
+            // faces only but include empty patches.
+
+                //- Numbering of boundary faces is face-mesh.nInternalFaces()
+                const globalIndex& globalBoundaryFaceNumbering() const;
+                const labelListList& globalPointBoundaryFaces() const;
+                const mapDistribute& globalPointBoundaryFacesMap() const;
+
+            // Coupled point to boundary cell
+
+                //- From boundary cell to mesh cell
+                const labelList& boundaryCells() const;
+
+                //- Numbering of boundary cells is according to boundaryCells()
+                const globalIndex& globalBoundaryCellNumbering() const;
+                const labelListList& globalPointBoundaryCells() const;
+                const mapDistribute& globalPointBoundaryCellsMap() const;
+
+
         // Edit
 
             //- Update for moving points.
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C
index bd94a30142..6dfc702fdd 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C
@@ -70,11 +70,11 @@ void Foam::globalPoints::addToSend
 (
     const primitivePatch& pp,
     const label patchPointI,
-    const procPointList& knownInfo,
+    const labelList& knownInfo,
 
     DynamicList<label>& patchFaces,
     DynamicList<label>& indexInFace,
-    DynamicList<procPointList>& allInfo
+    DynamicList<labelList>& allInfo
 )
 {
     label meshPointI = pp.meshPoints()[patchPointI];
@@ -97,58 +97,65 @@ void Foam::globalPoints::addToSend
 
 
 // Add nbrInfo to myInfo. Return true if anything changed.
-// nbrInfo is for a point a list of all the processors using it (and the
-// meshPoint label on that processor)
+// nbrInfo is for a point a list of all the global points using it
 bool Foam::globalPoints::mergeInfo
 (
-    const procPointList& nbrInfo,
-    procPointList& myInfo
+    const labelList& nbrInfo,
+    labelList& myInfo
 )
 {
-    // Indices of entries in nbrInfo not yet in myInfo.
-    DynamicList<label> newInfo(nbrInfo.size());
+    labelList newInfo(myInfo);
+    label newI = newInfo.size();
+    newInfo.setSize(newI + nbrInfo.size());
 
     forAll(nbrInfo, i)
     {
-        const procPoint& info = nbrInfo[i];
-
-        // Check if info already in myInfo.
-        label index = -1;
-
-        forAll(myInfo, j)
-        {
-            if (myInfo[j] == info)
-            {
-                // Already have information for processor/point combination
-                // in my list so skip.
-                index = j;
-
-                break;
-            }
-        }
+        label index = findIndex(myInfo, nbrInfo[i]);
 
         if (index == -1)
         {
-            // Mark this information as being not yet in myInfo
-            newInfo.append(i);
+            newInfo[newI++] = nbrInfo[i];
         }
     }
 
-    newInfo.shrink();
-
-    // Append all nbrInfos referenced in newInfo to myInfo.
-
-    label index = myInfo.size();
-
-    myInfo.setSize(index + newInfo.size());
-
-    forAll(newInfo, i)
-    {
-        myInfo[index++] = nbrInfo[newInfo[i]];
-    }
+    newInfo.setSize(newI);
 
     // Did anything change?
-    return newInfo.size() > 0;
+    bool anyChanged = (newI > myInfo.size());
+
+    myInfo.transfer(newInfo);
+
+    return anyChanged;
+}
+
+
+Foam::label Foam::globalPoints::meshToLocalPoint
+(
+    const Map<label>& meshToPatchPoint, // from mesh point to local numbering
+    const label meshPointI
+)
+{
+    return
+    (
+        meshToPatchPoint.size() == 0
+      ? meshPointI
+      : meshToPatchPoint[meshPointI]
+    );
+}
+
+
+Foam::label Foam::globalPoints::localToMeshPoint
+(
+    const labelList& patchToMeshPoint,
+    const label localPointI
+)
+{
+    return
+    (
+        patchToMeshPoint.size() == 0
+      ? localPointI
+      : patchToMeshPoint[localPointI]
+    );
 }
 
 
@@ -156,34 +163,30 @@ bool Foam::globalPoints::mergeInfo
 // Uses mergeInfo above. Returns true if data kept for meshPointI changed.
 bool Foam::globalPoints::storeInfo
 (
-    const procPointList& nbrInfo,
-    const label meshPointI
+    const labelList& nbrInfo,
+    const label localPointI
 )
 {
     label infoChanged = false;
 
     // Get the index into the procPoints list.
-    Map<label>::iterator iter = meshToProcPoint_.find(meshPointI);
+    Map<label>::iterator iter = meshToProcPoint_.find(localPointI);
 
     if (iter != meshToProcPoint_.end())
     {
-        procPointList& knownInfo = procPoints_[iter()];
-
-        if (mergeInfo(nbrInfo, knownInfo))
+        if (mergeInfo(nbrInfo, procPoints_[iter()]))
         {
             infoChanged = true;
         }
     }
     else
     {
-        procPointList knownInfo(1);
-        knownInfo[0][0] = Pstream::myProcNo();
-        knownInfo[0][1] = meshPointI;
+        labelList knownInfo(1, globalIndices_.toGlobal(localPointI));
 
         if (mergeInfo(nbrInfo, knownInfo))
         {
             // Update addressing from into procPoints
-            meshToProcPoint_.insert(meshPointI, procPoints_.size());
+            meshToProcPoint_.insert(localPointI, procPoints_.size());
             // Insert into list of equivalences.
             procPoints_.append(knownInfo);
 
@@ -197,6 +200,7 @@ bool Foam::globalPoints::storeInfo
 // Insert my own points into structure and mark as changed.
 void Foam::globalPoints::initOwnPoints
 (
+    const Map<label>& meshToPatchPoint,
     const bool allPoints,
     labelHashSet& changedPoints
 )
@@ -221,18 +225,24 @@ void Foam::globalPoints::initOwnPoints
                 forAll(meshPoints, i)
                 {
                     label meshPointI = meshPoints[i];
+                    label localPointI = meshToLocalPoint
+                    (
+                        meshToPatchPoint,
+                        meshPointI
+                    );
+                    labelList knownInfo
+                    (
+                        1,
+                        globalIndices_.toGlobal(localPointI)
+                    );
 
-                    procPointList knownInfo(1);
-                    knownInfo[0][0] = Pstream::myProcNo();
-                    knownInfo[0][1] = meshPointI;
-
-                    // Update addressing from meshpoint to index in procPoints
-                    meshToProcPoint_.insert(meshPointI, procPoints_.size());
+                    // Update addressing from point to index in procPoints
+                    meshToProcPoint_.insert(localPointI, procPoints_.size());
                     // Insert into list of equivalences.
                     procPoints_.append(knownInfo);
 
                     // Update changedpoints info.
-                    changedPoints.insert(meshPointI);
+                    changedPoints.insert(localPointI);
                 }
             }
             else
@@ -243,18 +253,25 @@ void Foam::globalPoints::initOwnPoints
                 forAll(boundaryPoints, i)
                 {
                     label meshPointI = meshPoints[boundaryPoints[i]];
+                    label localPointI = meshToLocalPoint
+                    (
+                        meshToPatchPoint,
+                        meshPointI
+                    );
 
-                    procPointList knownInfo(1);
-                    knownInfo[0][0] = Pstream::myProcNo();
-                    knownInfo[0][1] = meshPointI;
+                    labelList knownInfo
+                    (
+                        1,
+                        globalIndices_.toGlobal(localPointI)
+                    );
 
-                    // Update addressing from meshpoint to index in procPoints
-                    meshToProcPoint_.insert(meshPointI, procPoints_.size());
+                    // Update addressing from point to index in procPoints
+                    meshToProcPoint_.insert(localPointI, procPoints_.size());
                     // Insert into list of equivalences.
                     procPoints_.append(knownInfo);
 
                     // Update changedpoints info.
-                    changedPoints.insert(meshPointI);
+                    changedPoints.insert(localPointI);
                 }
             }
         }
@@ -263,8 +280,12 @@ void Foam::globalPoints::initOwnPoints
 
 
 // Send all my info on changedPoints_ to my neighbours.
-void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
- const
+void Foam::globalPoints::sendPatchPoints
+(
+    const Map<label>& meshToPatchPoint,
+    PstreamBuffers& pBufs,
+    const labelHashSet& changedPoints
+) const
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
@@ -280,10 +301,10 @@ void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
             // index in patch face
             DynamicList<label> indexInFace(pp.nPoints());
             // all information I currently hold about this patchPoint
-            DynamicList<procPointList> allInfo(pp.nPoints());
+            DynamicList<labelList> allInfo(pp.nPoints());
 
 
-            // Now collect information on all mesh points mentioned in
+            // Now collect information on all points mentioned in
             // changedPoints. Note that these points only should occur on
             // processorPatches (or rather this is a limitation!).
 
@@ -292,14 +313,19 @@ void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
             forAll(meshPoints, patchPointI)
             {
                 label meshPointI = meshPoints[patchPointI];
+                label localPointI = meshToLocalPoint
+                (
+                    meshToPatchPoint,
+                    meshPointI
+                );
 
-                if (changedPoints.found(meshPointI))
+                if (changedPoints.found(localPointI))
                 {
-                    label index = meshToProcPoint_[meshPointI];
+                    label index = meshToProcPoint_[localPointI];
 
-                    const procPointList& knownInfo = procPoints_[index];
+                    const labelList& knownInfo = procPoints_[index];
 
-                    // Add my information about meshPointI to the send buffers
+                    // Add my information about localPointI to the send buffers
                     addToSend
                     (
                         pp,
@@ -312,9 +338,6 @@ void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
                     );
                 }
             }
-            patchFaces.shrink();
-            indexInFace.shrink();
-            allInfo.shrink();
 
             // Send to neighbour
             {
@@ -328,12 +351,7 @@ void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
                         << patchFaces.size() << endl;
                 }
 
-                OPstream toNeighbour
-                (
-                    Pstream::blocking,
-                    procPatch.neighbProcNo()
-                );
-
+                UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
                 toNeighbour << patchFaces << indexInFace << allInfo;
             }
         }
@@ -345,8 +363,13 @@ void Foam::globalPoints::sendPatchPoints(const labelHashSet& changedPoints)
 // After finishing will have updated
 // - procPoints_ : all neighbour information merged in.
 // - meshToProcPoint_
-// - changedPoints: all meshPoints for which something changed.
-void Foam::globalPoints::receivePatchPoints(labelHashSet& changedPoints)
+// - changedPoints: all points for which something changed.
+void Foam::globalPoints::receivePatchPoints
+(
+    const Map<label>& meshToPatchPoint,
+    PstreamBuffers& pBufs,
+    labelHashSet& changedPoints
+)
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
@@ -364,14 +387,10 @@ void Foam::globalPoints::receivePatchPoints(labelHashSet& changedPoints)
 
             labelList patchFaces;
             labelList indexInFace;
-            List<procPointList> nbrInfo;
+            List<labelList> nbrInfo;
 
             {
-                IPstream fromNeighbour
-                (
-                    Pstream::blocking,
-                    procPatch.neighbProcNo()
-                );
+                UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
                 fromNeighbour >> patchFaces >> indexInFace >> nbrInfo;
             }
 
@@ -392,19 +411,15 @@ void Foam::globalPoints::receivePatchPoints(labelHashSet& changedPoints)
                 // Get the meshpoint on my side
                 label meshPointI = f[index];
 
-                //const procPointList& info = nbrInfo[i];
-                //Pout << "Received for my coord "
-                //    << mesh_.points()[meshPointI];
-                //
-                //forAll(info, j)
-                //{
-                //    Pout<< ' ' <<info[j];
-                //}
-                //Pout<< endl;
+                label localPointI = meshToLocalPoint
+                (
+                    meshToPatchPoint,
+                    meshPointI
+                );
 
-                if (storeInfo(nbrInfo[i], meshPointI))
+                if (storeInfo(nbrInfo[i], localPointI))
                 {
-                    changedPoints.insert(meshPointI);
+                    changedPoints.insert(localPointI);
                 }
             }
         }
@@ -428,29 +443,41 @@ void Foam::globalPoints::receivePatchPoints(labelHashSet& changedPoints)
                 label meshPointA = meshPoints[e[0]];
                 label meshPointB = meshPoints[e[1]];
 
-                // Do we have information on meshPointA?
+                label localA = meshToLocalPoint
+                (
+                    meshToPatchPoint,
+                    meshPointA
+                );
+                label localB = meshToLocalPoint
+                (
+                    meshToPatchPoint,
+                    meshPointB
+                );
+
+
+                // Do we have information on pointA?
                 Map<label>::iterator procPointA =
-                    meshToProcPoint_.find(meshPointA);
+                    meshToProcPoint_.find(localA);
 
                 if (procPointA != meshToProcPoint_.end())
                 {
                     // Store A info onto pointB
-                    if (storeInfo(procPoints_[procPointA()], meshPointB))
+                    if (storeInfo(procPoints_[procPointA()], localB))
                     {
-                        changedPoints.insert(meshPointB);
+                        changedPoints.insert(localB);
                     }
                 }
 
                 // Same for info on pointB
                 Map<label>::iterator procPointB =
-                    meshToProcPoint_.find(meshPointB);
+                    meshToProcPoint_.find(localB);
 
                 if (procPointB != meshToProcPoint_.end())
                 {
                     // Store B info onto pointA
-                    if (storeInfo(procPoints_[procPointB()], meshPointA))
+                    if (storeInfo(procPoints_[procPointB()], localA))
                     {
-                        changedPoints.insert(meshPointA);
+                        changedPoints.insert(localA);
                     }
                 }
             }
@@ -461,25 +488,24 @@ void Foam::globalPoints::receivePatchPoints(labelHashSet& changedPoints)
 
 // Remove entries which are handled by normal face-face communication. I.e.
 // those points where the equivalence list is only me and my (face)neighbour
-void Foam::globalPoints::remove(const Map<label>& directNeighbours)
+void Foam::globalPoints::remove
+(
+    const labelList& patchToMeshPoint,
+    const Map<label>& directNeighbours
+)
 {
     // Save old ones.
     Map<label> oldMeshToProcPoint(meshToProcPoint_);
     meshToProcPoint_.clear();
 
-    List<procPointList> oldProcPoints;
+    List<labelList> oldProcPoints;
     oldProcPoints.transfer(procPoints_);
 
     // Go through all equivalences
-    for
-    (
-        Map<label>::const_iterator iter = oldMeshToProcPoint.begin();
-        iter != oldMeshToProcPoint.end();
-        ++iter
-    )
+    forAllConstIter(Map<label>, oldMeshToProcPoint, iter)
     {
-        label meshPointI = iter.key();
-        const procPointList& pointInfo = oldProcPoints[iter()];
+        label localPointI = iter.key();
+        const labelList& pointInfo = oldProcPoints[iter()];
 
         if (pointInfo.size() == 2)
         {
@@ -488,23 +514,29 @@ void Foam::globalPoints::remove(const Map<label>& directNeighbours)
             // is in it. This would be an ordinary connection and can be
             // handled by normal face-face connectivity.
 
-            const procPoint& a = pointInfo[0];
-            const procPoint& b = pointInfo[1];
+            const label a = pointInfo[0];
+            const label b = pointInfo[1];
 
             if
             (
-                (a[0] == Pstream::myProcNo() && directNeighbours.found(a[1]))
-             || (b[0] == Pstream::myProcNo() && directNeighbours.found(b[1]))
+                (
+                    globalIndices_.isLocal(a)
+                 && directNeighbours.found(globalIndices_.toLocal(a))
+                )
+             || (
+                    globalIndices_.isLocal(b)
+                 && directNeighbours.found(globalIndices_.toLocal(b))
+                )
             )
             {
                 // Normal faceNeighbours
-                if (a[0] == Pstream::myProcNo())
+                if (globalIndices_.isLocal(a))
                 {
                     //Pout<< "Removing direct neighbour:"
                     //    << mesh_.points()[a[1]]
                     //    << endl;
                 }
-                else if (b[0] == Pstream::myProcNo())
+                else if (globalIndices_.isLocal(b))
                 {
                     //Pout<< "Removing direct neighbour:"
                     //    << mesh_.points()[b[1]]
@@ -525,7 +557,7 @@ void Foam::globalPoints::remove(const Map<label>& directNeighbours)
                 // be found if the two domains are face connected at all
                 // (not shown in the picture)
 
-                meshToProcPoint_.insert(meshPointI, procPoints_.size());
+                meshToProcPoint_.insert(localPointI, procPoints_.size());
                 procPoints_.append(pointInfo);
             }
         }
@@ -536,17 +568,17 @@ void Foam::globalPoints::remove(const Map<label>& directNeighbours)
             // So this meshPoint will have info of size one only.
             if
             (
-                pointInfo[0][0] != Pstream::myProcNo()
-             || !directNeighbours.found(pointInfo[0][1])
+                !globalIndices_.isLocal(pointInfo[0])
+             || !directNeighbours.found(globalIndices_.toLocal(pointInfo[0]))
             )
             {
-                meshToProcPoint_.insert(meshPointI, procPoints_.size());
+                meshToProcPoint_.insert(localPointI, procPoints_.size());
                 procPoints_.append(pointInfo);
             }
         }
         else
         {
-            meshToProcPoint_.insert(meshPointI, procPoints_.size());
+            meshToProcPoint_.insert(localPointI, procPoints_.size());
             procPoints_.append(pointInfo);
         }
     }
@@ -555,35 +587,67 @@ void Foam::globalPoints::remove(const Map<label>& directNeighbours)
 }
 
 
+// Compact indices
+void Foam::globalPoints::compact()
+{
+    // TBD: find same procPoints entries. Or rather check if
+    // in a procPoints there are points with indices < my index.
+    // This will just merge identical entries so lower storage, but will
+    // not affect anything else. Note: only relevant with cyclics.
+
+    labelList oldToNew(procPoints_.size(), -1);
+    labelList newToOld(meshToProcPoint_.size());
+
+    label newIndex = 0;
+    forAllIter(Map<label>, meshToProcPoint_, iter)
+    {
+        label oldIndex = iter();
+
+        if (oldToNew[oldIndex] == -1)
+        {
+            iter() = newIndex;
+            oldToNew[oldIndex] = newIndex;
+            newToOld[newIndex] = oldIndex;
+            newIndex++;
+        }
+    }
+    List<labelList> oldProcPoints;
+    oldProcPoints.transfer(procPoints_);
+
+    procPoints_.setSize(meshToProcPoint_.size());
+    forAll(procPoints_, i)
+    {
+        // Transfer
+        procPoints_[i].transfer(oldProcPoints[newToOld[i]]);
+    }
+}
+
+
 // Get (indices of) points for which I am master (= lowest numbered point on
 // lowest numbered processor).
 // (equivalence lists should be complete by now)
-Foam::labelList Foam::globalPoints::getMasterPoints() const
+Foam::labelList Foam::globalPoints::getMasterPoints
+(
+    const labelList& patchToMeshPoint
+) const
 {
     labelList masterPoints(nPatchPoints_);
     label nMaster = 0;
 
-    // Go through all equivalences and determine meshPoints where I am master.
-    for
-    (
-        Map<label>::const_iterator iter = meshToProcPoint_.begin();
-        iter != meshToProcPoint_.end();
-        ++iter
-    )
+    // Go through all equivalences and determine points where I am master.
+    forAllConstIter(Map<label>, meshToProcPoint_, iter)
     {
-        label meshPointI = iter.key();
-        const procPointList& pointInfo = procPoints_[iter()];
+        label localPointI = iter.key();
+        const labelList& pointInfo = procPoints_[iter()];
 
         if (pointInfo.size() < 2)
         {
             // Points should have an equivalence list >= 2 since otherwise
             // they would be direct neighbours and have been removed in the
             // call to 'remove'.
-            Pout<< "MeshPoint:" << meshPointI
-                << " coord:" << mesh_.points()[meshPointI]
-                << " has no corresponding point on a neighbouring processor"
-                << endl;
-            FatalErrorIn("globalPoints::getMasterPoints()")
+            label meshPointI = localToMeshPoint(patchToMeshPoint, localPointI);
+
+            FatalErrorIn("globalPoints::getMasterPoints(..)")
                 << '[' << Pstream::myProcNo() << ']'
                 << " MeshPoint:" << meshPointI
                 << " coord:" << mesh_.points()[meshPointI]
@@ -592,33 +656,17 @@ Foam::labelList Foam::globalPoints::getMasterPoints() const
         }
         else
         {
-            // Find lowest processor and lowest mesh point on this processor.
-            label lowestProcI = labelMax;
-            label lowestPointI = labelMax;
-
-            forAll(pointInfo, i)
-            {
-                label proc = pointInfo[i][0];
-
-                if (proc < lowestProcI)
-                {
-                    lowestProcI = proc;
-                    lowestPointI = pointInfo[i][1];
-                }
-                else if (proc == lowestProcI)
-                {
-                    lowestPointI = min(lowestPointI, pointInfo[i][1]);
-                }
-            }
+            // Check if lowest numbered processor and point is
+            // on this processor. Since already sorted is first element.
 
             if
             (
-                lowestProcI == Pstream::myProcNo()
-             && lowestPointI == meshPointI
+                globalIndices_.isLocal(pointInfo[0])
+             && globalIndices_.toLocal(pointInfo[0]) == localPointI
             )
             {
                 // I am lowest numbered processor and point. Add to my list.
-                masterPoints[nMaster++] = meshPointI;
+                masterPoints[nMaster++] = localPointI;
             }
         }
     }
@@ -630,7 +678,11 @@ Foam::labelList Foam::globalPoints::getMasterPoints() const
 
 
 // Send subset of lists
-void Foam::globalPoints::sendSharedPoints(const labelList& changedIndices) const
+void Foam::globalPoints::sendSharedPoints
+(
+    PstreamBuffers& pBufs,
+    const labelList& changedIndices
+) const
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
@@ -643,7 +695,7 @@ void Foam::globalPoints::sendSharedPoints(const labelList& changedIndices) const
             const processorPolyPatch& procPatch =
                 refCast<const processorPolyPatch>(pp);
 
-            OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
+            UOPstream toNeighbour(procPatch.neighbProcNo(), pBufs);
 
             if (debug)
             {
@@ -652,6 +704,7 @@ void Foam::globalPoints::sendSharedPoints(const labelList& changedIndices) const
                     << changedIndices.size() << endl;
             }
 
+            // Send over changed elements
             toNeighbour
                 << UIndirectList<label>(sharedPointAddr_, changedIndices)()
                 << UIndirectList<label>(sharedPointLabels_, changedIndices)();
@@ -661,10 +714,15 @@ void Foam::globalPoints::sendSharedPoints(const labelList& changedIndices) const
 
 
 // Receive shared point indices for all my shared points. Note that since
-// there are only a few here we can build a reverse map using the meshpoint
+// there are only a few here we can build a reverse map using the point label
 // instead of doing all this relative point indexing (patch face + index in
 // face) as in send/receivePatchPoints
-void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
+void Foam::globalPoints::receiveSharedPoints
+(
+    const Map<label>& meshToPatchPoint,
+    PstreamBuffers& pBufs,
+    labelList& changedIndices
+)
 {
     changedIndices.setSize(sharedPointAddr_.size());
     label nChanged = 0;
@@ -681,22 +739,18 @@ void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
             const processorPolyPatch& procPatch =
                 refCast<const processorPolyPatch>(pp);
 
-            // Map from neighbouring meshPoint to sharedPoint)
+            // Map from neighbouring mesh or patch point to sharedPoint)
             Map<label> nbrSharedPoints(sharedPointAddr_.size());
 
             {
-                // Receive meshPoints on neighbour and sharedPoints and build
+                // Receive points on neighbour and sharedPoints and build
                 // map from it. Note that we could have built the map on the
                 // neighbour and sent it over.
                 labelList nbrSharedPointAddr;
                 labelList nbrSharedPointLabels;
 
                 {
-                    IPstream fromNeighbour
-                    (
-                        Pstream::blocking,
-                        procPatch.neighbProcNo()
-                    );
+                    UIPstream fromNeighbour(procPatch.neighbProcNo(), pBufs);
                     fromNeighbour >> nbrSharedPointAddr >> nbrSharedPointLabels;
                 }
 
@@ -705,7 +759,7 @@ void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
                 {
                     nbrSharedPoints.insert
                     (
-                        nbrSharedPointLabels[i], // meshpoint on neighbour
+                        nbrSharedPointLabels[i], // mesh/patchpoint on neighbour
                         nbrSharedPointAddr[i]    // sharedPoint label
                     );
                 }
@@ -713,39 +767,36 @@ void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
 
 
             // Merge into whatever information I hold.
-            for
-            (
-                Map<label>::const_iterator iter = meshToProcPoint_.begin();
-                iter != meshToProcPoint_.end();
-                ++iter
-            )
+            forAllConstIter(Map<label>, meshToProcPoint_, iter)
             {
-                label meshPointI = iter.key();
+                label localPointI = iter.key();
                 label index = iter();
 
                 if (sharedPointAddr_[index] == -1)
                 {
-                    // No shared point known yet for this meshPoint.
+                    // No shared point known yet for this point.
                     // See if was received from neighbour.
-                    const procPointList& knownInfo = procPoints_[index];
+                    const labelList& knownInfo = procPoints_[index];
 
                     // Check through the whole equivalence list for any
                     // point from the neighbour.
                     forAll(knownInfo, j)
                     {
-                        const procPoint& info = knownInfo[j];
+                        const label info = knownInfo[j];
+                        label procI = globalIndices_.whichProcID(info);
+                        label pointI = globalIndices_.toLocal(procI, info);
 
                         if
                         (
-                            (info[0] == procPatch.neighbProcNo())
-                         && nbrSharedPoints.found(info[1])
+                            procI == procPatch.neighbProcNo()
+                         && nbrSharedPoints.found(pointI)
                         )
                         {
                             // So this knownInfo contains the neighbour point
-                            label sharedPointI = nbrSharedPoints[info[1]];
+                            label sharedPointI = nbrSharedPoints[pointI];
 
                             sharedPointAddr_[index] = sharedPointI;
-                            sharedPointLabels_[index] = meshPointI;
+                            sharedPointLabels_[index] = localPointI;
                             changedIndices[nChanged++] = index;
 
                             break;
@@ -759,13 +810,15 @@ void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
             const cyclicPolyPatch& cycPatch =
                 refCast<const cyclicPolyPatch>(pp);
 
-            // Build map from meshPoint to sharedPoint
-            Map<label> meshToSharedPoint(sharedPointAddr_.size());
+            // Build map from mesh or patch point to sharedPoint
+            Map<label> localToSharedPoint(sharedPointAddr_.size());
             forAll(sharedPointLabels_, i)
             {
-                label meshPointI = sharedPointLabels_[i];
-
-                meshToSharedPoint.insert(meshPointI, sharedPointAddr_[i]);
+                localToSharedPoint.insert
+                (
+                    sharedPointLabels_[i],
+                    sharedPointAddr_[i]
+                );
             }
 
             // Sync all info.
@@ -775,17 +828,27 @@ void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
             forAll(connections, i)
             {
                 const edge& e = connections[i];
-
                 label meshPointA = pp.meshPoints()[e[0]];
                 label meshPointB = pp.meshPoints()[e[1]];
 
-                // Do we already have shared point for meshPointA?
-                Map<label>::iterator fndA = meshToSharedPoint.find(meshPointA);
-                Map<label>::iterator fndB = meshToSharedPoint.find(meshPointB);
+                label localA = meshToLocalPoint
+                (
+                    meshToPatchPoint,
+                    meshPointA
+                );
+                label localB = meshToLocalPoint
+                (
+                    meshToPatchPoint,
+                    meshPointB
+                );
 
-                if (fndA != meshToSharedPoint.end())
+                // Do we already have shared point for pointA?
+                Map<label>::iterator fndA = localToSharedPoint.find(localA);
+                Map<label>::iterator fndB = localToSharedPoint.find(localB);
+
+                if (fndA != localToSharedPoint.end())
                 {
-                    if (fndB != meshToSharedPoint.end())
+                    if (fndB != localToSharedPoint.end())
                     {
                         if (fndA() != fndB())
                         {
@@ -808,26 +871,26 @@ void Foam::globalPoints::receiveSharedPoints(labelList& changedIndices)
                         // No shared point yet for B.
                         label sharedPointI = fndA();
 
-                        // Store shared point for meshPointB
-                        label index = meshToProcPoint_[meshPointB];
+                        // Store shared point for pointB
+                        label index = meshToProcPoint_[localB];
 
                         sharedPointAddr_[index] = sharedPointI;
-                        sharedPointLabels_[index] = meshPointB;
+                        sharedPointLabels_[index] = localB;
                         changedIndices[nChanged++] = index;
                     }
                 }
                 else
                 {
                     // No shared point yet for A.
-                    if (fndB != meshToSharedPoint.end())
+                    if (fndB != localToSharedPoint.end())
                     {
                         label sharedPointI = fndB();
 
-                        // Store shared point for meshPointA
-                        label index = meshToProcPoint_[meshPointA];
+                        // Store shared point for pointA
+                        label index = meshToProcPoint_[localA];
 
                         sharedPointAddr_[index] = sharedPointI;
-                        sharedPointLabels_[index] = meshPointA;
+                        sharedPointLabels_[index] = localA;
                         changedIndices[nChanged++] = index;
                     }
                 }
@@ -887,18 +950,12 @@ Foam::edgeList Foam::globalPoints::coupledPoints(const cyclicPolyPatch& pp)
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-Foam::globalPoints::globalPoints(const polyMesh& mesh)
-:
-    mesh_(mesh),
-    nPatchPoints_(countPatchPoints(mesh.boundaryMesh())),
-    procPoints_(nPatchPoints_),
-    meshToProcPoint_(nPatchPoints_),
-    sharedPointAddr_(0),
-    sharedPointLabels_(0),
-    nGlobalPoints_(0)
+void Foam::globalPoints::calculateSharedPoints
+(
+    const Map<label>& meshToPatchPoint, // from mesh point to local numbering
+    const labelList& patchToMeshPoint,  // from local numbering to mesh point
+    const bool keepAllPoints
+)
 {
     if (debug)
     {
@@ -921,24 +978,33 @@ Foam::globalPoints::globalPoints(const polyMesh& mesh)
     //   shared point is not on the boundary of any processor patches using it.
     //   This would happen if a domain was pinched such that two patches share
     //   a point or edge.
-    initOwnPoints(true, changedPoints);
+    initOwnPoints(meshToPatchPoint, true, changedPoints);
 
     // Do one exchange iteration to get neighbour points.
-    sendPatchPoints(changedPoints);
-    receivePatchPoints(changedPoints);
+    {
+        PstreamBuffers pBufs(Pstream::defaultCommsType);
+        sendPatchPoints(meshToPatchPoint, pBufs, changedPoints);
+        pBufs.finishedSends();
+        receivePatchPoints(meshToPatchPoint, pBufs, changedPoints);
+    }
 
 
     // Save neighbours reachable through face-face communication.
-    Map<label> neighbourList(meshToProcPoint_);
-
+    Map<label> neighbourList;
+    if (!keepAllPoints)
+    {
+        neighbourList = meshToProcPoint_;
+    }
 
     // Exchange until nothing changes on all processors.
     bool changed = false;
 
     do
     {
-        sendPatchPoints(changedPoints);
-        receivePatchPoints(changedPoints);
+        PstreamBuffers pBufs(Pstream::defaultCommsType);
+        sendPatchPoints(meshToPatchPoint, pBufs, changedPoints);
+        pBufs.finishedSends();
+        receivePatchPoints(meshToPatchPoint, pBufs, changedPoints);
 
         changed = changedPoints.size() > 0;
         reduce(changed, orOp<bool>());
@@ -947,29 +1013,57 @@ Foam::globalPoints::globalPoints(const polyMesh& mesh)
 
 
     // Remove direct neighbours from point equivalences.
-    remove(neighbourList);
+    if (!keepAllPoints)
+    {
+        remove(patchToMeshPoint, neighbourList);
+    }
+    else
+    {
+        // Compact out unused elements
+        compact();
+    }
+
+    procPoints_.shrink();
+
+    // Sort procPoints in incremental order. This will make the master the
+    // first element.
+    forAllIter(Map<label>, meshToProcPoint_, iter)
+    {
+        sort(procPoints_[iter()]);
+    }
 
 
-    //Pout<< "After removing locally connected points:" << endl;
-    //for
-    //(
-    //    Map<label>::const_iterator iter = meshToProcPoint_.begin();
-    //    iter != meshToProcPoint_.end();
-    //    ++iter
-    //)
-    //{
-    //    label meshPointI = iter.key();
-    //    const procPointList& pointInfo = procPoints_[iter()];
-    //
-    //    forAll(pointInfo, i)
-    //    {
-    //        Pout<< "    pointI:" << meshPointI << ' '
-    //            << mesh.points()[meshPointI]
-    //            << " connected to proc " << pointInfo[i][0]
-    //            << " point:" << pointInfo[i][1]
-    //        << endl;
-    //    }
-    //}
+    // Pout<< "Now connected points:" << endl;
+    // forAllConstIter(Map<label>, meshToProcPoint_, iter)
+    // {
+    //     label localI = iter.key();
+    //     const labelList& pointInfo = procPoints_[iter()];
+    // 
+    //     Pout<< "pointI:" << localI << " index:" << iter()
+    //         << " coord:"
+    //         << mesh_.points()[localToMeshPoint(patchToMeshPoint, localI)]
+    //         << endl;
+    // 
+    //     forAll(pointInfo, i)
+    //     {
+    //         label procI = globalIndices_.whichProcID(pointInfo[i]);
+    //         Pout<< "    connected to proc " << procI
+    //             << " localpoint:"
+    //             << globalIndices_.toLocal(procI, pointInfo[i]);
+    // 
+    //         if (globalIndices_.isLocal(pointInfo[i]))
+    //         {
+    //             label meshPointI = localToMeshPoint
+    //             (
+    //                 patchToMeshPoint,
+    //                 globalIndices_.toLocal(pointInfo[i])
+    //             );
+    //             Pout<< " at:" <<  mesh_.points()[meshPointI];
+    //         }
+    // 
+    //         Pout<< endl;
+    //     }
+    // }
 
 
     // We now have - in procPoints_ - a list of points which are shared between
@@ -988,49 +1082,21 @@ Foam::globalPoints::globalPoints(const polyMesh& mesh)
     sharedPointLabels_ = -1;
 
 
-    // Get points for which I am master (lowest numbered proc)
-    labelList masterPoints(getMasterPoints());
+    // Get point labels of points for which I am master (lowest
+    // numbered proc)
+    labelList masterPoints(getMasterPoints(patchToMeshPoint));
 
-
-    // Determine number of master points on all processors.
-    labelList sharedPointSizes(Pstream::nProcs());
-    sharedPointSizes[Pstream::myProcNo()] = masterPoints.size();
-
-    Pstream::gatherList(sharedPointSizes);
-    Pstream::scatterList(sharedPointSizes);
-
-    if (debug)
-    {
-        Pout<< "sharedPointSizes:" << sharedPointSizes << endl;
-    }
-
-    // Get total number of shared points
-    nGlobalPoints_ = 0;
-    forAll(sharedPointSizes, procI)
-    {
-        nGlobalPoints_ += sharedPointSizes[procI];
-    }
-
-    // Assign sharedPoint labels. Every processor gets assigned consecutive
-    // numbers for its master points.
-    // These are assigned in processor order so processor0 gets
-    // 0..masterPoints.size()-1 etc.
-
-    // My point labels start after those of lower numbered processors
-    label sharedPointI = 0;
-    for (label procI = 0; procI < Pstream::myProcNo(); procI++)
-    {
-        sharedPointI += sharedPointSizes[procI];
-    }
+    // Get global numbering for master points
+    globalIndex globalMasterPoints(masterPoints.size());
+    nGlobalPoints_ = globalMasterPoints.size();
 
     forAll(masterPoints, i)
     {
-        label meshPointI = masterPoints[i];
+        label localPointI = masterPoints[i];
+        label index = meshToProcPoint_[localPointI];
 
-        label index = meshToProcPoint_[meshPointI];
-
-        sharedPointLabels_[index] = meshPointI;
-        sharedPointAddr_[index] = sharedPointI++;
+        sharedPointLabels_[index] = localPointI;
+        sharedPointAddr_[index] = globalMasterPoints.toGlobal(i);
     }
 
 
@@ -1060,8 +1126,10 @@ Foam::globalPoints::globalPoints(const polyMesh& mesh)
             Pout<< "Determined " << changedIndices.size() << " shared points."
                 << " Exchanging them" << endl;
         }
-        sendSharedPoints(changedIndices);
-        receiveSharedPoints(changedIndices);
+        PstreamBuffers pBufs(Pstream::defaultCommsType);
+        sendSharedPoints(pBufs, changedIndices);
+        pBufs.finishedSends();
+        receiveSharedPoints(meshToPatchPoint, pBufs, changedIndices);
 
         changed = changedIndices.size() > 0;
         reduce(changed, orOp<bool>());
@@ -1073,7 +1141,7 @@ Foam::globalPoints::globalPoints(const polyMesh& mesh)
     {
         if (sharedPointLabels_[i] == -1)
         {
-            FatalErrorIn("globalPoints::globalPoints(const polyMesh& mesh)")
+            FatalErrorIn("globalPoints::globalPoints(const polyMesh&)")
                 << "Problem: shared point on processor " << Pstream::myProcNo()
                 << " not set at index " << sharedPointLabels_[i] << endl
                 << "This might mean the individual processor domains are not"
@@ -1091,4 +1159,56 @@ Foam::globalPoints::globalPoints(const polyMesh& mesh)
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from mesh
+Foam::globalPoints::globalPoints
+(
+    const polyMesh& mesh,
+    const bool keepAllPoints
+)
+:
+    mesh_(mesh),
+    globalIndices_(mesh_.nPoints()),
+    nPatchPoints_(countPatchPoints(mesh.boundaryMesh())),
+    procPoints_(nPatchPoints_),
+    meshToProcPoint_(nPatchPoints_),
+    sharedPointAddr_(0),
+    sharedPointLabels_(0),
+    nGlobalPoints_(0)
+{
+    // Empty patch maps to signal storing mesh point labels
+    Map<label> meshToPatchPoint(0);
+    labelList patchToMeshPoint(0);
+
+    calculateSharedPoints(meshToPatchPoint, patchToMeshPoint, keepAllPoints);
+}
+
+
+// Construct from mesh and patch of coupled faces
+Foam::globalPoints::globalPoints
+(
+    const polyMesh& mesh,
+    const indirectPrimitivePatch& coupledPatch,
+    const bool keepAllPoints
+)
+:
+    mesh_(mesh),
+    globalIndices_(coupledPatch.nPoints()),
+    nPatchPoints_(coupledPatch.nPoints()),
+    procPoints_(nPatchPoints_),
+    meshToProcPoint_(nPatchPoints_),
+    sharedPointAddr_(0),
+    sharedPointLabels_(0),
+    nGlobalPoints_(0)
+{
+    calculateSharedPoints
+    (
+        coupledPatch.meshPointMap(),
+        coupledPatch.meshPoints(),
+        keepAllPoints
+    );
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H
index c2361da2ba..41d010ec9f 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.H
@@ -29,9 +29,8 @@ Description
     Calculates points shared by more than two processor patches or cyclic
     patches.
 
-    Is used in globalMeshData. (this info is needed for point-edge
-    communication where you do all but these shared points with patch to
-    patch communication but need to do a reduce on these shared points)
+    Is used in globalMeshData. (this info is needed for point/edge
+    communication where processor swaps are not enough to exchange data)
 
     Works purely topological and using local communication only.
     Needs:
@@ -41,8 +40,8 @@ Description
       - f[0] ordering on patch faces to be ok.
 
     Works by constructing equivalence lists for all the points on processor
-    patches. These list are procPointList and give processor and meshPoint
-    label on that processor.
+    patches. These list are in globalIndex numbering (so consecutively numbered
+    per processor)
     E.g.
     @verbatim
           ((7 93)(4 731)(3 114))
@@ -85,6 +84,9 @@ Description
      endloop until nothing changes.
 
 
+    Note: the data held is either mesh point labels (construct from mesh only)
+    or patch point labels (construct from mesh and patch).
+
 SourceFiles
     globalPoints.C
 
@@ -95,11 +97,11 @@ SourceFiles
 
 #include "DynamicList.H"
 #include "Map.H"
-#include "labelList.H"
-#include "FixedList.H"
 #include "primitivePatch.H"
 #include "className.H"
 #include "edgeList.H"
+#include "globalIndex.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -119,30 +121,29 @@ class globalPoints
 {
     // Private classes
 
-        //- Define procPointList as holding a list of meshPoint/processor labels
-        typedef FixedList<label, 2> procPoint;
-        typedef List<procPoint> procPointList;
-
     // Private data
 
         //- Mesh reference
         const polyMesh& mesh_;
 
+        //- Global numbering of points
+        globalIndex globalIndices_;
+
         //- Sum of points on processor patches (unfiltered, point on 2 patches
         //  counts as 2)
         const label nPatchPoints_;
 
         //- All points on boundaries and their corresponding connected points
         //  on other processors.
-        DynamicList<procPointList> procPoints_;
+        DynamicList<labelList> procPoints_;
 
-        //- Map from mesh point to index in procPoints
+        //- Map from mesh (or patch) point to index in procPoints
         Map<label> meshToProcPoint_;
 
         //- Shared points used by this processor (= global point number)
         labelList sharedPointAddr_;
 
-        //- My meshpoints corresponding to the shared points
+        //- My mesh(or patch) points corresponding to the shared points
         labelList sharedPointLabels_;
 
         //- Total number of shared points.
@@ -155,56 +156,108 @@ class globalPoints
         //  information is collected.
         static label countPatchPoints(const polyBoundaryMesh&);
 
+        ////- Get all faces on coupled patches
+        //static labelListl coupledFaces(const polyBoundaryMesh&);
+
         //- Add information about patchPointI in relative indices to send
         //  buffers (patchFaces, indexInFace etc.)
         static void addToSend
         (
             const primitivePatch&,
             const label patchPointI,
-            const procPointList&,
+            const labelList&,
             DynamicList<label>& patchFaces,
             DynamicList<label>& indexInFace,
-            DynamicList<procPointList>& allInfo
+            DynamicList<labelList>& allInfo
         );
 
         //- Merge info from neighbour into my data
         static bool mergeInfo
         (
-            const procPointList& nbrInfo,
-            procPointList& myInfo
+            const labelList& nbrInfo,
+            labelList& myInfo
+        );
+
+        //- From mesh point to 'local point'. Is the mesh point itself
+        //  if meshToPatchPoint is empty.
+        static label meshToLocalPoint
+        (
+            const Map<label>& meshToPatchPoint,
+            const label meshPointI
+        );
+
+        //- Opposite of meshToLocalPoint.
+        static label localToMeshPoint
+        (
+            const labelList& patchToMeshPoint,
+            const label localPointI
         );
 
         //- Store (and merge) info for meshPointI
-        bool storeInfo(const procPointList& nbrInfo, const label meshPointI);
+        bool storeInfo
+        (
+            const labelList& nbrInfo,
+            const label localPointI
+        );
 
         //- Initialize procPoints_ to my patch points. allPoints = true:
         //  seed with all patch points, = false: only boundaryPoints().
-        void initOwnPoints(const bool allPoints, labelHashSet& changedPoints);
+        void initOwnPoints
+        (
+            const Map<label>& meshToPatchPoint,
+            const bool allPoints,
+            labelHashSet& changedPoints
+        );
 
         //- Send subset of procPoints to neighbours
-        void sendPatchPoints(const labelHashSet& changedPoints) const;
+        void sendPatchPoints
+        (
+            const Map<label>&,
+            PstreamBuffers&,
+            const labelHashSet&
+        ) const;
 
         //- Receive neighbour points and merge into my procPoints.
-        void receivePatchPoints(labelHashSet& changedPoints);
+        void receivePatchPoints
+        (
+            const Map<label>&,
+            PstreamBuffers&,
+            labelHashSet&
+        );
 
         //- Remove entries of size 2 where meshPoint is in provided Map.
         //  Used to remove normal face-face connected points.
-        void remove(const Map<label>&);
+        void remove(const labelList& patchToMeshPoint, const Map<label>&);
+
+        //- Compact out unused elements of procPoints.
+        void compact();
 
         //- Get indices of point for which I am master (lowest numbered proc)
-        labelList getMasterPoints() const;
+        labelList getMasterPoints(const labelList& patchToMeshPoint) const;
 
         //- Send subset of shared points to neighbours
-        void sendSharedPoints(const labelList& changedIndices) const;
+        void sendSharedPoints(PstreamBuffers&, const labelList&) const;
 
         //- Receive shared points and update subset.
-        void receiveSharedPoints(labelList& changedIndices);
+        void receiveSharedPoints
+        (
+            const Map<label>&,
+            PstreamBuffers&,
+            labelList&
+        );
 
-
-        //- Should move into cyclicPolyPatch but some ordering problem
+        //- Should move into cyclicPolyPatch ordering problem
         //  keeps on giving problems.
         static edgeList coupledPoints(const cyclicPolyPatch&);
 
+        //- Do all calculations.
+        void calculateSharedPoints
+        (
+            const Map<label>&,
+            const labelList&,
+            const bool keepAllPoints
+        );
+
         //- Disallow default bitwise copy construct
         globalPoints(const globalPoints&);
 
@@ -220,36 +273,59 @@ public:
 
     // Constructors
 
-        //- Construct from mesh
-        globalPoints(const polyMesh& mesh);
+        //- Construct from mesh.
+        //  keepAllPoints = false : filter out points that are on two
+        //  neighbouring coupled patches (so can be swapped)
+        globalPoints(const polyMesh& mesh, const bool keepAllPoints);
+
+        //- Construct from mesh and patch of coupled faces. Difference with
+        //  construct from mesh only is that this stores the meshToProcPoint,
+        //  procPoints and sharedPointLabels as patch local point labels
+        //  instead of mesh point labels.
+        globalPoints
+        (
+            const polyMesh& mesh,
+            const indirectPrimitivePatch& coupledPatch,
+            const bool keepAllPoints
+        );
 
 
     // Member Functions
 
         // Access
 
-            label nPatchPoints() const
-            {
-                return nPatchPoints_;
-            }
-
+            //- From (mesh or patch) point to index in procPoints
             const Map<label>& meshToProcPoint() const
             {
                 return meshToProcPoint_;
             }
 
+            //- procPoints is per point the connected points (in global
+            //  point numbers)
+            const DynamicList<labelList>& procPoints() const
+            {
+                return procPoints_;
+            }
+
+            //- Global numbering of (mesh or patch) points
+            const globalIndex& globalIndices() const
+            {
+                return globalIndices_;
+            }
+
             //- shared points used by this processor (= global point number)
             const labelList& sharedPointAddr() const
             {
                 return sharedPointAddr_;
             }
 
-            //- my meshpoints corresponding to the shared points
+            //- my (mesh or patch)points corresponding to the shared points
             const labelList& sharedPointLabels() const
             {
                 return sharedPointLabels_;
             }
 
+            //- total number of shared points
             label nGlobalPoints() const
             {
                 return nGlobalPoints_;
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C
index 47620111b8..f9b4a37e46 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C
@@ -288,6 +288,7 @@ Foam::mapDistribute::mapDistribute
 
         forAll(compactMap, procI)
         {
+            compactMap[procI].clear();
             if (procI != Pstream::myProcNo())
             {
                 compactMap[procI].resize(2*nNonLocal[procI]);
@@ -433,6 +434,7 @@ Foam::mapDistribute::mapDistribute
 
         forAll(compactMap, procI)
         {
+            compactMap[procI].clear();
             if (procI != Pstream::myProcNo())
             {
                 compactMap[procI].resize(2*nNonLocal[procI]);
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
index e0af4d6a0c..b1f78236de 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.H
@@ -276,6 +276,49 @@ public:
                 }
             }
 
+            //- Reverse distribute data using default commsType.
+            template<class T>
+            void reverseDistribute(const label constructSize, List<T>& fld)
+            const
+            {
+                if (Pstream::defaultCommsType == Pstream::nonBlocking)
+                {
+                    distribute
+                    (
+                        Pstream::nonBlocking,
+                        List<labelPair>(),
+                        constructSize,
+                        constructMap_,
+                        subMap_,
+                        fld
+                    );
+                }
+                else if (Pstream::defaultCommsType == Pstream::scheduled)
+                {
+                    distribute
+                    (
+                        Pstream::scheduled,
+                        schedule(),
+                        constructSize,
+                        constructMap_,
+                        subMap_,
+                        fld
+                    );
+                }
+                else
+                {
+                    distribute
+                    (
+                        Pstream::blocking,
+                        List<labelPair>(),
+                        constructSize,
+                        constructMap_,
+                        subMap_,
+                        fld
+                    );
+                }
+            }
+
             //- Correct for topo change.
             void updateMesh(const mapPolyMesh&)
             {
diff --git a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
index 09ce0cb1b7..29d9fc5c96 100644
--- a/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C
@@ -552,7 +552,6 @@ void Foam::mapDistribute::distribute
     {
         if (!contiguous<T>())
         {
-//XXXXXX
             PstreamBuffers pBuffs(Pstream::nonBlocking);
 
             // Stream data into buffer
diff --git a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
index 754eaf907c..ddd0b18cb2 100644
--- a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
+++ b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
@@ -26,7 +26,7 @@ Class
     Foam::syncTools
 
 Description
-    Various tools to aid synchronizing lists across coupled patches.
+    Various tools to aid synchronizing lists across coupled patches. WIP.
 
     Require
     - combineOperator (e.g. sumEqOp - not sumOp!) that is defined for the
@@ -34,8 +34,11 @@ Description
     - null value which gets overridden by any valid value.
     - transform function
 
-    Can apply coordinate rotation/separation on cyclics but only for faces
+    note:Can apply coordinate rotation/separation on cyclics but only for faces
     or if there is a single rotation/separation tensor.
+    note:syncPointList or syncEdgeList will visit shared points/edges
+    multiple times (once through patch exchange, once through shared
+    points reduce). Should be replaced by pointMesh functionality.
 
 SourceFiles
     syncTools.C
diff --git a/src/OpenFOAM/primitives/chars/wchar/wchar.H b/src/OpenFOAM/primitives/chars/wchar/wchar.H
index c731a3dc81..61519a51d3 100644
--- a/src/OpenFOAM/primitives/chars/wchar/wchar.H
+++ b/src/OpenFOAM/primitives/chars/wchar/wchar.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,10 @@ Description
 SourceFiles
     wcharIO.C
 
+SeeAlso
+    http://en.wikipedia.org/wiki/UTF-8
+    http://en.wikibooks.org/wiki/Unicode/Character_reference
+
 \*---------------------------------------------------------------------------*/
 
 #ifndef wchar_H
@@ -48,10 +52,10 @@ class Ostream;
 
 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
 
-//- Output via a normal char
+//- Output wide character (Unicode) as UTF-8
 Ostream& operator<<(Ostream&, const wchar_t);
 
-//- Output string via normal char
+//- Output wide character (Unicode) string as UTF-8
 Ostream& operator<<(Ostream&, const wchar_t*);
 
 
diff --git a/src/OpenFOAM/primitives/chars/wchar/wcharIO.C b/src/OpenFOAM/primitives/chars/wchar/wcharIO.C
index 47a5ed0560..9796d4b6eb 100644
--- a/src/OpenFOAM/primitives/chars/wchar/wcharIO.C
+++ b/src/OpenFOAM/primitives/chars/wchar/wcharIO.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -22,11 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-    Reads a char from an input stream, for a given version
-    number and File format. If an ascii File is being read, then the line
-    numbers are counted and an erroneous read is reported.
-
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
@@ -38,7 +33,68 @@ Description
 
 Foam::Ostream& Foam::operator<<(Ostream& os, const wchar_t wc)
 {
-    os.write(char(wc));
+    if (!(wc & ~0x0000007F))
+    {
+        // 0x00000000 - 0x0000007F: (1-byte output)
+        // 0xxxxxxx
+        os.write(char(wc));
+    }
+    else if (!(wc & ~0x000007FF))
+    {
+        // 0x00000080 - 0x000007FF: (2-byte output)
+        // 110bbbaa 10aaaaaa
+        os.write(char(0xC0 | ((wc >> 6) & 0x1F)));
+        os.write(char(0x80 | ((wc) & 0x3F)));
+    }
+    else if (!(wc & ~0x0000FFFF))
+    {
+        // 0x00000800 - 0x0000FFFF: (3-byte output)
+        // 1110bbbb 10bbbbaa 10aaaaaa
+        os.write(char(0xE0 | ((wc >> 12) & 0x0F)));
+        os.write(char(0x80 | ((wc >> 6) & 0x3F)));
+        os.write(char(0x80 | ((wc) & 0x3F)));
+    }
+    else if (!(wc & ~0x001FFFFF))
+    {
+        // 0x00010000 - 0x001FFFFF: (4-byte output)
+        // 11110ccc 10ccbbbb 10bbbbaa 10aaaaaa
+        os.write(char(0xF0 | ((wc >> 18) & 0x07)));
+        os.write(char(0x80 | ((wc >> 12) & 0x3F)));
+        os.write(char(0x80 | ((wc >> 6) & 0x3F)));
+        os.write(char(0x80 | ((wc) & 0x3F)));
+    }
+    else if (!(wc & ~0x03FFFFFF))
+    {
+        // 0x00200000 - 0x03FFFFFF: (5-byte output)
+        // 111110dd 10cccccc 10ccbbbb 10bbbbaa 10aaaaaa
+        os.write(char(0xF8 | ((wc >> 24) & 0x03)));
+        os.write(char(0x80 | ((wc >> 18) & 0x3F)));
+        os.write(char(0x80 | ((wc >> 12) & 0x3F)));
+        os.write(char(0x80 | ((wc >> 6) & 0x3F)));
+        os.write(char(0x80 | ((wc) & 0x3F)));
+    }
+    else if (!(wc & ~0x7FFFFFFF))
+    {
+        // 0x04000000 - 0x7FFFFFFF: (6-byte output)
+        // 1111110d 10dddddd 10cccccc 10ccbbbb 10bbbbaa 10aaaaaa
+        os.write(char(0xFC | ((wc >> 30) & 0x01)));
+        os.write(char(0x80 | ((wc >> 24) & 0x3F)));
+        os.write(char(0x80 | ((wc >> 18) & 0x3F)));
+        os.write(char(0x80 | ((wc >> 12) & 0x3F)));
+        os.write(char(0x80 | ((wc >> 6) & 0x3F)));
+        os.write(char(0x80 | ((wc) & 0x3F)));
+    }
+    else
+    {
+        // according to man page utf8(7)
+        // the Unicode standard specifies no characters above 0x0010FFFF,
+        // so Unicode characters can only be up to four bytes long in UTF-8.
+
+        // report anything unknown/invalid as replacement character U+FFFD
+        os.write(char(0xEF));
+        os.write(char(0xBF));
+        os.write(char(0xBD));
+    }
 
     os.check("Ostream& operator<<(Ostream&, const wchar_t)");
     return os;
@@ -51,10 +107,10 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const wchar_t* ws)
     {
         for (const wchar_t* p = ws; *p; ++p)
         {
-            os.write(char(*p));
+            os  << *p;
         }
     }
-    os.check("Ostream& operator<<(Ostream&, const wchar_t*)");
+
     return os;
 }
 
diff --git a/src/edgeMesh/Make/options b/src/edgeMesh/Make/options
index e69de29bb2..7e207d0dbe 100644
--- a/src/edgeMesh/Make/options
+++ b/src/edgeMesh/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude
+
+LIB_LIBS = \
+    -lfileFormats
diff --git a/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C b/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
index 4b178d1dbd..fec0c883a3 100644
--- a/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
+++ b/src/edgeMesh/edgeFormats/nas/NASedgeFormat.C
@@ -29,34 +29,6 @@ License
 #include "IStringStream.H"
 #include "PackedBoolList.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-// Do weird things to extract a floating point number
-Foam::scalar Foam::fileFormats::NASedgeFormat::parseNASCoord
-(
-    const string& s
-)
-{
-    size_t expSign = s.find_last_of("+-");
-
-    if (expSign != string::npos && expSign > 0 && !isspace(s[expSign-1]))
-    {
-        scalar mantissa = readScalar(IStringStream(s.substr(0, expSign))());
-        scalar exponent = readScalar(IStringStream(s.substr(expSign+1))());
-
-        if (s[expSign] == '-')
-        {
-            exponent = -exponent;
-        }
-        return mantissa * pow(10, exponent);
-    }
-    else
-    {
-        return readScalar(IStringStream(s)());
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fileFormats::NASedgeFormat::NASedgeFormat
diff --git a/src/edgeMesh/edgeFormats/nas/NASedgeFormat.H b/src/edgeMesh/edgeFormats/nas/NASedgeFormat.H
index 600b880cb7..3ccd4397af 100644
--- a/src/edgeMesh/edgeFormats/nas/NASedgeFormat.H
+++ b/src/edgeMesh/edgeFormats/nas/NASedgeFormat.H
@@ -37,6 +37,7 @@ SourceFiles
 #define NASedgeFormat_H
 
 #include "edgeMesh.H"
+#include "NASCore.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,7 +52,8 @@ namespace fileFormats
 
 class NASedgeFormat
 :
-    public edgeMesh
+    public edgeMesh,
+    public NASCore
 {
     // Private Member Functions
 
@@ -61,13 +63,6 @@ class NASedgeFormat
         //- Disallow default bitwise assignment
         void operator=(const NASedgeFormat&);
 
-protected:
-
-    // Protected Member Functions
-
-        //- Do weird things to extract number
-        static scalar parseNASCoord(const string&);
-
 public:
 
     // Constructors
diff --git a/src/edgeMesh/edgeFormats/starcd/STARCDedgeFormat.C b/src/edgeMesh/edgeFormats/starcd/STARCDedgeFormat.C
index b6ecd4d19e..f2fec25cc5 100644
--- a/src/edgeMesh/edgeFormats/starcd/STARCDedgeFormat.C
+++ b/src/edgeMesh/edgeFormats/starcd/STARCDedgeFormat.C
@@ -59,134 +59,6 @@ inline void Foam::fileFormats::STARCDedgeFormat::writeLines
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bool Foam::fileFormats::STARCDedgeFormat::readHeader
-(
-    IFstream& is,
-    const word& signature
-)
-{
-    if (!is.good())
-    {
-        FatalErrorIn
-        (
-            "fileFormats::STARCDedgeFormat::readHeader(...)"
-        )
-            << "cannot read " << signature  << "  " << is.name()
-            << abort(FatalError);
-    }
-
-    word header;
-    label majorVersion;
-
-    string line;
-
-    is.getLine(line);
-    IStringStream(line)() >> header;
-
-    is.getLine(line);
-    IStringStream(line)() >> majorVersion;
-
-    // add other checks ...
-    if (header != signature)
-    {
-        Info<< "header mismatch " << signature << "  " << is.name()
-            << endl;
-    }
-
-    return true;
-}
-
-
-void Foam::fileFormats::STARCDedgeFormat::writeHeader
-(
-    Ostream& os,
-    const char* filetype
-)
-{
-    os  << "PROSTAR_" << filetype << nl
-        << 4000
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << endl;
-}
-
-
-bool Foam::fileFormats::STARCDedgeFormat::readPoints
-(
-    IFstream& is,
-    pointField& points,
-    labelList& ids
-)
-{
-    //
-    // read .vrt file
-    // ~~~~~~~~~~~~~~
-
-    if (!is.good())
-    {
-        FatalErrorIn
-        (
-            "fileFormats::STARCDedgeFormat::readPoints(...)"
-        )
-            << "Cannot read file " << is.name()
-            << exit(FatalError);
-    }
-
-    readHeader(is, "PROSTAR_VERTEX");
-
-    DynamicList<point> dynPoints;
-    DynamicList<label> dynPointId;    // STAR-CD index of points
-
-    label lineLabel;
-    while ((is >> lineLabel).good())
-    {
-        scalar x, y, z;
-
-        is >> x >> y >> z;
-
-        dynPoints.append(point(x, y, z));
-        dynPointId.append(lineLabel);
-    }
-
-    points.transfer(dynPoints);
-    ids.transfer(dynPointId);
-
-    return true;
-}
-
-
-
-void Foam::fileFormats::STARCDedgeFormat::writePoints
-(
-    Ostream& os,
-    const pointField& pointLst
-)
-{
-    writeHeader(os, "VERTEX");
-
-    // Set the precision of the points data to 10
-    os.precision(10);
-
-    // force decimal point for Fortran input
-    os.setf(std::ios::showpoint);
-
-    forAll(pointLst, ptI)
-    {
-        os
-            << ptI + 1 << " "
-            << pointLst[ptI].x() << " "
-            << pointLst[ptI].y() << " "
-            << pointLst[ptI].z() << nl;
-    }
-    os.flush();
-}
-
-
 void Foam::fileFormats::STARCDedgeFormat::writeCase
 (
     Ostream& os,
diff --git a/src/edgeMesh/edgeFormats/starcd/STARCDedgeFormat.H b/src/edgeMesh/edgeFormats/starcd/STARCDedgeFormat.H
index 31d4a8db5d..d8ffb6b790 100644
--- a/src/edgeMesh/edgeFormats/starcd/STARCDedgeFormat.H
+++ b/src/edgeMesh/edgeFormats/starcd/STARCDedgeFormat.H
@@ -43,6 +43,8 @@ SourceFiles
 #define STARCDedgeFormat_H
 
 #include "edgeMesh.H"
+#include "STARCDCore.H"
+
 #include "IFstream.H"
 #include "Ostream.H"
 #include "OFstream.H"
@@ -60,7 +62,8 @@ namespace fileFormats
 
 class STARCDedgeFormat
 :
-    public edgeMesh
+    public edgeMesh,
+    public STARCDCore
 {
     // Private Data
 
@@ -90,14 +93,6 @@ protected:
 
     // Protected Member Functions
 
-    static bool readHeader(IFstream&, const word&);
-
-    static void writeHeader(Ostream&, const char* filetype);
-
-    static bool readPoints(IFstream&, pointField&, labelList& ids);
-
-    static void writePoints(Ostream&, const pointField&);
-
     static void writeCase
     (
         Ostream&,
diff --git a/src/fileFormats/Make/files b/src/fileFormats/Make/files
new file mode 100644
index 0000000000..a3acab834d
--- /dev/null
+++ b/src/fileFormats/Make/files
@@ -0,0 +1,4 @@
+nas/NASCore.C
+starcd/STARCDCore.C
+
+LIB = $(FOAM_LIBBIN)/libfileFormats
diff --git a/src/fileFormats/Make/options b/src/fileFormats/Make/options
new file mode 100644
index 0000000000..e69de29bb2
diff --git a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormatCore.C b/src/fileFormats/nas/NASCore.C
similarity index 87%
rename from src/surfMesh/surfaceFormats/nas/NASsurfaceFormatCore.C
rename to src/fileFormats/nas/NASCore.C
index 97b6a9ae40..9af082eef1 100644
--- a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormatCore.C
+++ b/src/fileFormats/nas/NASCore.C
@@ -24,13 +24,18 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "NASsurfaceFormatCore.H"
+#include "NASCore.H"
 #include "IStringStream.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Do weird things to extract a floating point number
-Foam::scalar Foam::fileFormats::NASsurfaceFormatCore::parseNASCoord
+Foam::fileFormats::NASCore::NASCore()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::fileFormats::NASCore::parseNASCoord
 (
     const string& s
 )
diff --git a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormatCore.H b/src/fileFormats/nas/NASCore.H
similarity index 81%
rename from src/surfMesh/surfaceFormats/nas/NASsurfaceFormatCore.H
rename to src/fileFormats/nas/NASCore.H
index 2c7b32d781..cbafd16de9 100644
--- a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormatCore.H
+++ b/src/fileFormats/nas/NASCore.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,21 +23,21 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::fileFormats::NASsurfaceFormatCore
+    Foam::fileFormats::NASCore
 
 Description
-    Internal class used by the NASsurfaceFormat
+    Core routines used when reading/writing NASTRAN files.
 
 SourceFiles
-    NASsurfaceFormatCore.C
+    NASCore.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef NASsurfaceFormatCore_H
-#define NASsurfaceFormatCore_H
+#ifndef NASCore_H
+#define NASCore_H
 
-#include "Ostream.H"
-#include "OFstream.H"
+#include "scalar.H"
+#include "string.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,17 +47,25 @@ namespace fileFormats
 {
 
 /*---------------------------------------------------------------------------*\
-                    Class NASsurfaceFormatCore Declaration
+                    Class fileFormats::NASCore Declaration
 \*---------------------------------------------------------------------------*/
 
-class NASsurfaceFormatCore
+class NASCore
 {
-protected:
+public:
 
-    // Protected Member Functions
+    // Public Member Functions
 
         //- Do weird things to extract number
         static scalar parseNASCoord(const string&);
+
+
+    // Constructors
+
+        //- Construct null
+        NASCore();
+
+
 };
 
 
diff --git a/src/fileFormats/starcd/STARCDCore.C b/src/fileFormats/starcd/STARCDCore.C
new file mode 100644
index 0000000000..e8206c29d8
--- /dev/null
+++ b/src/fileFormats/starcd/STARCDCore.C
@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "STARCDCore.H"
+#include "ListOps.H"
+#include "clock.H"
+#include "PackedBoolList.H"
+#include "IStringStream.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fileFormats::STARCDCore::STARCDCore()
+{}
+
+
+// * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
+
+bool Foam::fileFormats::STARCDCore::readHeader
+(
+    IFstream& is,
+    const word& signature
+)
+{
+    if (!is.good())
+    {
+        FatalErrorIn
+        (
+            "fileFormats::STARCDCore::readHeader(...)"
+        )
+            << "cannot read " << signature  << "  " << is.name()
+            << abort(FatalError);
+    }
+
+    word header;
+    label majorVersion;
+
+    string line;
+
+    is.getLine(line);
+    IStringStream(line)() >> header;
+
+    is.getLine(line);
+    IStringStream(line)() >> majorVersion;
+
+    // add other checks ...
+    if (header != signature)
+    {
+        Info<< "header mismatch " << signature << "  " << is.name()
+            << endl;
+    }
+
+    return true;
+}
+
+
+void Foam::fileFormats::STARCDCore::writeHeader
+(
+    Ostream& os,
+    const word& filetype
+)
+{
+    os  << "PROSTAR_" << filetype << nl
+        << 4000
+        << " " << 0
+        << " " << 0
+        << " " << 0
+        << " " << 0
+        << " " << 0
+        << " " << 0
+        << " " << 0
+        << endl;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::fileFormats::STARCDCore::readPoints
+(
+    IFstream& is,
+    pointField& points,
+    labelList& ids
+)
+{
+    if (!is.good())
+    {
+        FatalErrorIn
+        (
+            "fileFormats::STARCDedgeFormat::readPoints(...)"
+        )
+            << "Cannot read file " << is.name()
+            << exit(FatalError);
+    }
+
+    readHeader(is, "PROSTAR_VERTEX");
+
+
+    // reuse memory if possible
+    DynamicList<point> dynPoints(points.xfer());
+    DynamicList<label> dynPointId(ids.xfer());    // STAR-CD index of points
+
+    dynPoints.clear();
+    dynPointId.clear();
+
+    label lineLabel;
+    while ((is >> lineLabel).good())
+    {
+        scalar x, y, z;
+
+        is >> x >> y >> z;
+
+        dynPoints.append(point(x, y, z));
+        dynPointId.append(lineLabel);
+    }
+
+    points.transfer(dynPoints);
+    ids.transfer(dynPointId);
+
+    return true;
+}
+
+
+void Foam::fileFormats::STARCDCore::writePoints
+(
+    Ostream& os,
+    const pointField& pointLst
+)
+{
+    writeHeader(os, "VERTEX");
+
+    // Set the precision of the points data to 10
+    os.precision(10);
+
+    // force decimal point for Fortran input
+    os.setf(std::ios::showpoint);
+
+    forAll(pointLst, ptI)
+    {
+        os
+            << ptI + 1 << " "
+            << pointLst[ptI].x() << " "
+            << pointLst[ptI].y() << " "
+            << pointLst[ptI].z() << nl;
+    }
+    os.flush();
+}
+
+
+
+
+// ************************************************************************* //
diff --git a/src/fileFormats/starcd/STARCDCore.H b/src/fileFormats/starcd/STARCDCore.H
new file mode 100644
index 0000000000..7057d842cb
--- /dev/null
+++ b/src/fileFormats/starcd/STARCDCore.H
@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2006-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::fileFormats::STARCDCore
+
+Description
+    Core routines used when reading/writing pro-STAR vrt/cel/bnd files.
+
+SourceFiles
+    STARCDCore.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef STARCDCore_H
+#define STARCDCore_H
+
+#include "IFstream.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace fileFormats
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class fileFormats::STARCDCore Declaration
+\*---------------------------------------------------------------------------*/
+
+class STARCDCore
+{
+protected:
+
+    // Protected Member Functions
+
+        //- Read header
+        static bool readHeader(IFstream&, const word& fileSignature);
+
+        //- Write header for fileType (CELL|VERTEX|BOUNDARY)
+        static void writeHeader(Ostream&, const word& fileType);
+
+
+protected:
+
+        enum cellType
+        {
+            starcdFluidType   = 1,
+            starcdSolidType   = 2,
+            starcdBaffleType  = 3,
+            starcdShellType   = 4,
+            starcdLineType    = 5,
+            starcdPointType   = 6
+        };
+
+        enum shapeType
+        {
+            starcdPoint = 1,
+            starcdLine  = 2,
+            starcdShell = 3,
+            starcdHex   = 11,
+            starcdPrism = 12,
+            starcdTet   = 13,
+            starcdPyr   = 14,
+            starcdPoly  = 255
+        };
+
+
+public:
+
+    // Public Member Functions
+
+        //- Read points from a (.vrt) file
+        // The file format is as follows:
+        // @verbatim
+        // Line 1:
+        //   PROSTAR_VERTEX  newline
+        //
+        // Line 2:
+        //   {version} 0 0 0 0 0 0 0  newline
+        //
+        // Body:
+        //   {vertexId}  {x}  {y}  {z}  newline
+        // @endverbatim
+        static bool readPoints(IFstream&, pointField&, labelList& ids);
+
+        //- Write header and points to (.vrt) file
+        static void writePoints(Ostream&, const pointField&);
+
+
+    // Constructors
+
+        //- Construct null
+        STARCDCore();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fileFormats
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C
index 8643dba429..ef24cbf2aa 100644
--- a/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/basic/sliced/slicedFvPatchField.C
@@ -55,7 +55,7 @@ slicedFvPatchField<Type>::slicedFvPatchField
     const DimensionedField<Type, volMesh>& iF
 )
 :
-    fvPatchField<Type>(p, iF)
+    fvPatchField<Type>(p, iF, Field<Type>())
 {}
 
 
diff --git a/src/mesh/autoMesh/Make/options b/src/mesh/autoMesh/Make/options
index 573d0007c1..28c3d5fdf0 100644
--- a/src/mesh/autoMesh/Make/options
+++ b/src/mesh/autoMesh/Make/options
@@ -8,7 +8,6 @@ EXE_INC = \
     -I$(LIB_SRC)/triSurface/lnInclude
 
 LIB_LIBS = \
-    -ldecompositionMethods \
     -ldynamicMesh \
     -lfiniteVolume \
     -llagrangian \
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
index 83852a1aff..5bfa725225 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoSnapDriver.C
@@ -1506,8 +1506,12 @@ void Foam::autoSnapDriver::doSnap
         // meshMover.
         calcNearestSurface(snapDist, meshMover);
 
-        // Get smoothly varying internal displacement field.
-        smoothDisplacement(snapParams, meshMover);
+        //// Get smoothly varying internal displacement field.
+        //- 2009-12-16 : was not found to be beneficial. Keeping internal
+        // fields fixed slightly increases skewness (on boundary)
+        // but lowers non-orthogonality quite a bit (e.g. 65->59 degrees).
+        // Maybe if better smoother?
+        //smoothDisplacement(snapParams, meshMover);
 
         // Apply internal displacement to mesh.
         scaleMesh(snapParams, nInitErrors, baffles, meshMover);
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index ada844b5a5..30e68c56a7 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -9,16 +9,17 @@ cellDist/wallPoint/wallPoint.C
 
 cellFeatures/cellFeatures.C
 
-coordinateSystems/coordinateSystem.C
-coordinateSystems/coordinateSystemNew.C
-coordinateSystems/coordinateSystems.C
-coordinateSystems/cylindricalCS.C
-coordinateSystems/sphericalCS.C
-coordinateSystems/parabolicCylindricalCS.C
-coordinateSystems/toroidalCS.C
-coordinateSystems/coordinateRotation/coordinateRotation.C
-coordinateSystems/coordinateRotation/EulerCoordinateRotation.C
-coordinateSystems/coordinateRotation/STARCDCoordinateRotation.C
+csys = coordinateSystems
+$(csys)/coordinateSystem.C
+$(csys)/coordinateSystemNew.C
+$(csys)/coordinateSystems.C
+$(csys)/cylindricalCS.C
+$(csys)/sphericalCS.C
+$(csys)/parabolicCylindricalCS.C
+$(csys)/toroidalCS.C
+$(csys)/coordinateRotation/coordinateRotation.C
+$(csys)/coordinateRotation/EulerCoordinateRotation.C
+$(csys)/coordinateRotation/STARCDCoordinateRotation.C
 
 edgeFaceCirculator/edgeFaceCirculator.C
 
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 0ff2202926..6bf989d4f5 100644
--- a/src/meshTools/Make/options
+++ b/src/meshTools/Make/options
@@ -6,4 +6,5 @@ EXE_INC = \
 LIB_LIBS = \
     -ltriSurface \
     -ldecompositionMethods \
+    -lmetisDecompositionMethod \
     -llagrangian
diff --git a/src/meshTools/indexedOctree/indexedOctree.H b/src/meshTools/indexedOctree/indexedOctree.H
index 6c2bfa3c6e..b653dc3e79 100644
--- a/src/meshTools/indexedOctree/indexedOctree.H
+++ b/src/meshTools/indexedOctree/indexedOctree.H
@@ -85,8 +85,7 @@ public:
         };
 
 
-    // Tree node. Has up pointer and down pointers.
-
+        //- Tree node. Has up pointer and down pointers.
         class node
         {
         public:
diff --git a/src/parallel/Allwmake b/src/parallel/Allwmake
index 34d323251a..48c7d41624 100755
--- a/src/parallel/Allwmake
+++ b/src/parallel/Allwmake
@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso decompositionMethods
+wmake libso metisDecomp
 wmake libso reconstruct
 
 if [ -d "$FOAM_MPI_LIBBIN" ]
diff --git a/src/parallel/decompositionMethods/Make/files b/src/parallel/decompositionMethods/Make/files
index 2acbfd5606..751ebe90f2 100644
--- a/src/parallel/decompositionMethods/Make/files
+++ b/src/parallel/decompositionMethods/Make/files
@@ -6,6 +6,4 @@ manualDecomp/manualDecomp.C
 
 scotchDecomp/scotchDecomp.C
 
-metisDecomp/metisDecomp.C
-
 LIB = $(FOAM_LIBBIN)/libdecompositionMethods
diff --git a/src/parallel/decompositionMethods/Make/options b/src/parallel/decompositionMethods/Make/options
index 97569aa663..fe6e8b6a05 100644
--- a/src/parallel/decompositionMethods/Make/options
+++ b/src/parallel/decompositionMethods/Make/options
@@ -1,8 +1,6 @@
 EXE_INC = \
-    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/src/libscotch/lnInclude \
-    -I$(WM_THIRD_PARTY_DIR)/metis-5.0pre2/include
+    -I$(WM_THIRD_PARTY_DIR)/scotch_5.1/include \
+    -I../metisDecomp/lnInclude
 
 LIB_LIBS = \
-    -lscotch \
-    -lmetis \
-    -lGKlib
+    -lscotch -lscotcherrexit
diff --git a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C
index 61a9483837..5375fb529b 100644
--- a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C
+++ b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.C
@@ -61,17 +61,12 @@ License
 #include "scotchDecomp.H"
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
-#include "IFstream.H"
 #include "Time.H"
 #include "cyclicPolyPatch.H"
 #include "OFstream.H"
-#include "metisDecomp.H"
 
 extern "C"
 {
-#define OMPI_SKIP_MPICXX
-#include "module.h"
-#include "common.h"
 #include "scotch.h"
 }
 
@@ -115,13 +110,13 @@ void Foam::scotchDecomp::check(const int retVal, const char* str)
     if (retVal)
     {
         FatalErrorIn("scotchDecomp::decompose(..)")
-            << "Called to scotch routine " << str << " failed."
+            << "Call to scotch routine " << str << " failed."
             << exit(FatalError);
     }
 }
 
 
-// Call Metis with options from dictionary.
+// Call scotch with options from dictionary.
 Foam::label Foam::scotchDecomp::decompose
 (
     const List<int>& adjncy,
@@ -173,7 +168,7 @@ Foam::label Foam::scotchDecomp::decompose
         {
             FatalErrorIn
             (
-                "parMetisDecomp::decompose"
+                "scotchDecomp::decompose"
                 "(const pointField&, const scalarField&)"
             )   << "Number of cell weights " << cWeights.size()
                 << " does not equal number of cells " << xadj.size()-1
@@ -377,7 +372,10 @@ Foam::labelList Foam::scotchDecomp::decompose
 {
     if (points.size() != mesh_.nCells())
     {
-        FatalErrorIn("scotchDecomp::decompose(const pointField&)")
+        FatalErrorIn
+        (
+            "scotchDecomp::decompose(const pointField&, const scalarField&)"
+        )
             << "Can use this decomposition method only for the whole mesh"
             << endl
             << "and supply one coordinate (cellCentre) for every cell." << endl
@@ -391,12 +389,7 @@ Foam::labelList Foam::scotchDecomp::decompose
     //   xadj(celli) : start of information in adjncy for celli
     List<int> adjncy;
     List<int> xadj;
-    metisDecomp::calcMetisCSR
-    (
-        mesh_,
-        adjncy,
-        xadj
-    );
+    calcCSR(mesh_, adjncy, xadj);
 
     // Decompose using default weights
     List<int> finalDecomp;
@@ -445,7 +438,7 @@ Foam::labelList Foam::scotchDecomp::decompose
             cellCells
         );
 
-        metisDecomp::calcMetisCSR(cellCells, adjncy, xadj);
+        calcCSR(cellCells, adjncy, xadj);
     }
 
     // Decompose using weights
@@ -475,7 +468,8 @@ Foam::labelList Foam::scotchDecomp::decompose
     {
         FatalErrorIn
         (
-            "scotchDecomp::decompose(const pointField&, const labelListList&)"
+            "scotchDecomp::decompose"
+            "(const labelListList&, const pointField&, const scalarField&)"
         )   << "Inconsistent number of cells (" << globalCellCells.size()
             << ") and number of cell centres (" << cellCentres.size()
             << ")." << exit(FatalError);
@@ -488,7 +482,7 @@ Foam::labelList Foam::scotchDecomp::decompose
 
     List<int> adjncy;
     List<int> xadj;
-    metisDecomp::calcMetisCSR(globalCellCells, adjncy, xadj);
+    calcCSR(globalCellCells, adjncy, xadj);
 
     // Decompose using weights
     List<int> finalDecomp;
@@ -504,4 +498,144 @@ Foam::labelList Foam::scotchDecomp::decompose
 }
 
 
+void Foam::scotchDecomp::calcCSR
+(
+    const polyMesh& mesh,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Make Metis CSR (Compressed Storage Format) storage
+    //   adjncy      : contains neighbours (= edges in graph)
+    //   xadj(celli) : start of information in adjncy for celli
+
+    xadj.setSize(mesh.nCells()+1);
+
+    // Initialise the number of internal faces of the cells to twice the
+    // number of internal faces
+    label nInternalFaces = 2*mesh.nInternalFaces();
+
+    // Check the boundary for coupled patches and add to the number of
+    // internal faces
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    forAll(pbm, patchi)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        {
+            nInternalFaces += pbm[patchi].size();
+        }
+    }
+
+    // Create the adjncy array the size of the total number of internal and
+    // coupled faces
+    adjncy.setSize(nInternalFaces);
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+    label freeAdj = 0;
+
+    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
+    {
+        xadj[cellI] = freeAdj;
+
+        const labelList& cFaces = mesh.cells()[cellI];
+
+        forAll(cFaces, i)
+        {
+            label faceI = cFaces[i];
+
+            if
+            (
+                mesh.isInternalFace(faceI)
+             || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
+            )
+            {
+                freeAdj++;
+            }
+        }
+    }
+    xadj[mesh.nCells()] = freeAdj;
+
+
+    // Fill in adjncy
+    // ~~~~~~~~~~~~~~
+
+    labelList nFacesPerCell(mesh.nCells(), 0);
+
+    // Internal faces
+    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+    {
+        label own = mesh.faceOwner()[faceI];
+        label nei = mesh.faceNeighbour()[faceI];
+
+        adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+    }
+
+    // Coupled faces. Only cyclics done.
+    forAll(pbm, patchi)
+    {
+        if (isA<cyclicPolyPatch>(pbm[patchi]))
+        {
+            const unallocLabelList& faceCells = pbm[patchi].faceCells();
+
+            label sizeby2 = faceCells.size()/2;
+
+            for (label facei=0; facei<sizeby2; facei++)
+            {
+                label own = faceCells[facei];
+                label nei = faceCells[facei + sizeby2];
+
+                adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
+                adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
+            }
+        }
+    }
+}
+
+
+// From cell-cell connections to Metis format (like CompactListList)
+void Foam::scotchDecomp::calcCSR
+(
+    const labelListList& cellCells,
+    List<int>& adjncy,
+    List<int>& xadj
+)
+{
+    // Count number of internal faces
+    label nConnections = 0;
+
+    forAll(cellCells, coarseI)
+    {
+        nConnections += cellCells[coarseI].size();
+    }
+
+    // Create the adjncy array as twice the size of the total number of
+    // internal faces
+    adjncy.setSize(nConnections);
+
+    xadj.setSize(cellCells.size()+1);
+
+
+    // Fill in xadj
+    // ~~~~~~~~~~~~
+    label freeAdj = 0;
+
+    forAll(cellCells, coarseI)
+    {
+        xadj[coarseI] = freeAdj;
+
+        const labelList& cCells = cellCells[coarseI];
+
+        forAll(cCells, i)
+        {
+            adjncy[freeAdj++] = cCells[i];
+        }
+    }
+    xadj[cellCells.size()] = freeAdj;
+}
+
+
+
 // ************************************************************************* //
diff --git a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H
index 54efb1bafa..93bcc7b9a4 100644
--- a/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H
+++ b/src/parallel/decompositionMethods/scotchDecomp/scotchDecomp.H
@@ -135,6 +135,25 @@ public:
             const scalarField& cWeights
         );
 
+
+        //- Helper to convert local connectivity (supplied as owner,neighbour)
+        //  into CSR (Metis,scotch) storage. Does cyclics but not processor
+        //  patches
+        static void calcCSR
+        (
+            const polyMesh& mesh,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
+
+        //- Helper to convert connectivity (supplied as cellcells) into
+        //  CSR (Metis,scotch) storage
+        static void calcCSR
+        (
+            const labelListList& globalCellCells,
+            List<int>& adjncy,
+            List<int>& xadj
+        );
 };
 
 
diff --git a/src/parallel/metisDecomp/Make/files b/src/parallel/metisDecomp/Make/files
new file mode 100644
index 0000000000..eb1e633e18
--- /dev/null
+++ b/src/parallel/metisDecomp/Make/files
@@ -0,0 +1,3 @@
+metisDecomp.C
+
+LIB = $(FOAM_LIBBIN)/libmetisDecompositionMethod
diff --git a/src/parallel/metisDecomp/Make/options b/src/parallel/metisDecomp/Make/options
new file mode 100644
index 0000000000..a7d5398f03
--- /dev/null
+++ b/src/parallel/metisDecomp/Make/options
@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(WM_THIRD_PARTY_DIR)/metis-5.0pre2/include \
+    -I../decompositionMethods/lnInclude
+
+LIB_LIBS = \
+    -lmetis \
+    -lGKlib
diff --git a/src/parallel/decompositionMethods/metisDecomp/metisDecomp.C b/src/parallel/metisDecomp/metisDecomp.C
similarity index 78%
rename from src/parallel/decompositionMethods/metisDecomp/metisDecomp.C
rename to src/parallel/metisDecomp/metisDecomp.C
index 3b2c5fb58c..9066e92248 100644
--- a/src/parallel/decompositionMethods/metisDecomp/metisDecomp.C
+++ b/src/parallel/metisDecomp/metisDecomp.C
@@ -27,9 +27,8 @@ License
 #include "metisDecomp.H"
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
-#include "IFstream.H"
 #include "Time.H"
-#include "cyclicPolyPatch.H"
+#include "scotchDecomp.H"
 
 extern "C"
 {
@@ -331,12 +330,7 @@ Foam::labelList Foam::metisDecomp::decompose
 
     List<int> adjncy;
     List<int> xadj;
-    calcMetisCSR
-    (
-        mesh_,
-        adjncy,
-        xadj
-    );
+    scotchDecomp::calcCSR(mesh_, adjncy, xadj);
 
     // Decompose using default weights
     List<int> finalDecomp;
@@ -352,145 +346,6 @@ Foam::labelList Foam::metisDecomp::decompose
 }
 
 
-void Foam::metisDecomp::calcMetisCSR
-(
-    const polyMesh& mesh,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    // Make Metis CSR (Compressed Storage Format) storage
-    //   adjncy      : contains neighbours (= edges in graph)
-    //   xadj(celli) : start of information in adjncy for celli
-
-    xadj.setSize(mesh.nCells()+1);
-
-    // Initialise the number of internal faces of the cells to twice the
-    // number of internal faces
-    label nInternalFaces = 2*mesh.nInternalFaces();
-
-    // Check the boundary for coupled patches and add to the number of
-    // internal faces
-    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
-
-    forAll(pbm, patchi)
-    {
-        if (isA<cyclicPolyPatch>(pbm[patchi]))
-        {
-            nInternalFaces += pbm[patchi].size();
-        }
-    }
-
-    // Create the adjncy array the size of the total number of internal and
-    // coupled faces
-    adjncy.setSize(nInternalFaces);
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-    label freeAdj = 0;
-
-    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
-    {
-        xadj[cellI] = freeAdj;
-
-        const labelList& cFaces = mesh.cells()[cellI];
-
-        forAll(cFaces, i)
-        {
-            label faceI = cFaces[i];
-
-            if
-            (
-                mesh.isInternalFace(faceI)
-             || isA<cyclicPolyPatch>(pbm[pbm.whichPatch(faceI)])
-            )
-            {
-                freeAdj++;
-            }
-        }
-    }
-    xadj[mesh.nCells()] = freeAdj;
-
-
-    // Fill in adjncy
-    // ~~~~~~~~~~~~~~
-
-    labelList nFacesPerCell(mesh.nCells(), 0);
-
-    // Internal faces
-    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-    {
-        label own = mesh.faceOwner()[faceI];
-        label nei = mesh.faceNeighbour()[faceI];
-
-        adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-        adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
-    }
-
-    // Coupled faces. Only cyclics done.
-    forAll(pbm, patchi)
-    {
-        if (isA<cyclicPolyPatch>(pbm[patchi]))
-        {
-            const unallocLabelList& faceCells = pbm[patchi].faceCells();
-
-            label sizeby2 = faceCells.size()/2;
-
-            for (label facei=0; facei<sizeby2; facei++)
-            {
-                label own = faceCells[facei];
-                label nei = faceCells[facei + sizeby2];
-
-                adjncy[xadj[own] + nFacesPerCell[own]++] = nei;
-                adjncy[xadj[nei] + nFacesPerCell[nei]++] = own;
-            }
-        }
-    }
-}
-
-
-// From cell-cell connections to Metis format (like CompactListList)
-void Foam::metisDecomp::calcMetisCSR
-(
-    const labelListList& cellCells,
-    List<int>& adjncy,
-    List<int>& xadj
-)
-{
-    // Count number of internal faces
-    label nConnections = 0;
-
-    forAll(cellCells, coarseI)
-    {
-        nConnections += cellCells[coarseI].size();
-    }
-
-    // Create the adjncy array as twice the size of the total number of
-    // internal faces
-    adjncy.setSize(nConnections);
-
-    xadj.setSize(cellCells.size()+1);
-
-
-    // Fill in xadj
-    // ~~~~~~~~~~~~
-    label freeAdj = 0;
-
-    forAll(cellCells, coarseI)
-    {
-        xadj[coarseI] = freeAdj;
-
-        const labelList& cCells = cellCells[coarseI];
-
-        forAll(cCells, i)
-        {
-            adjncy[freeAdj++] = cCells[i];
-        }
-    }
-    xadj[cellCells.size()] = freeAdj;
-}
-
-
 Foam::labelList Foam::metisDecomp::decompose
 (
     const labelList& agglom,
@@ -525,7 +380,7 @@ Foam::labelList Foam::metisDecomp::decompose
             cellCells
         );
 
-        calcMetisCSR(cellCells, adjncy, xadj);
+        scotchDecomp::calcCSR(cellCells, adjncy, xadj);
     }
 
     // Decompose using default weights
@@ -570,7 +425,7 @@ Foam::labelList Foam::metisDecomp::decompose
 
     List<int> adjncy;
     List<int> xadj;
-    calcMetisCSR(globalCellCells, adjncy, xadj);
+    scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
 
 
     // Decompose using default weights
diff --git a/src/parallel/decompositionMethods/metisDecomp/metisDecomp.H b/src/parallel/metisDecomp/metisDecomp.H
similarity index 89%
rename from src/parallel/decompositionMethods/metisDecomp/metisDecomp.H
rename to src/parallel/metisDecomp/metisDecomp.H
index 6a38e25c8a..c585686f60 100644
--- a/src/parallel/decompositionMethods/metisDecomp/metisDecomp.H
+++ b/src/parallel/metisDecomp/metisDecomp.H
@@ -132,24 +132,6 @@ public:
             const scalarField& cWeights
         );
 
-        //- Helper to convert connectivity (supplied as owner,neighbour) into
-        //  Metis storage
-        static void calcMetisCSR
-        (
-            const polyMesh& mesh,
-            List<int>& adjncy,
-            List<int>& xadj
-        );
-
-        //- Helper to convert connectivity (supplied as cellcells) into
-        //  Metis storage
-        static void calcMetisCSR
-        (
-            const labelListList& globalCellCells,
-            List<int>& adjncy,
-            List<int>& xadj
-        );
-
 };
 
 
diff --git a/src/parallel/parMetisDecomp/Make/options b/src/parallel/parMetisDecomp/Make/options
index 475bf5dbb1..79beb369af 100644
--- a/src/parallel/parMetisDecomp/Make/options
+++ b/src/parallel/parMetisDecomp/Make/options
@@ -4,7 +4,8 @@ EXE_INC = \
     $(PFLAGS) $(PINC) \
     -I$(WM_THIRD_PARTY_DIR)/ParMetis-3.1/ParMETISLib \
     -I$(WM_THIRD_PARTY_DIR)/ParMetis-3.1 \
-    -I../decompositionMethods/lnInclude
+    -I../decompositionMethods/lnInclude \
+    -I../metisDecomp/lnInclude
 
 LIB_LIBS = \
     -L$(FOAM_MPI_LIBBIN) \
diff --git a/src/parallel/parMetisDecomp/parMetisDecomp.C b/src/parallel/parMetisDecomp/parMetisDecomp.C
index c4dde21d65..cf7434b5c9 100644
--- a/src/parallel/parMetisDecomp/parMetisDecomp.C
+++ b/src/parallel/parMetisDecomp/parMetisDecomp.C
@@ -26,6 +26,7 @@ License
 
 #include "parMetisDecomp.H"
 #include "metisDecomp.H"
+#include "scotchDecomp.H"
 #include "syncTools.H"
 #include "addToRunTimeSelectionTable.H"
 #include "floatScalar.H"
@@ -434,7 +435,7 @@ Foam::labelList Foam::parMetisDecomp::decompose
         {
             FatalErrorIn
             (
-                "metisDecomp::decompose"
+                "parMetisDecomp::decompose"
                 "(const pointField&, const scalarField&)"
             )   << "Number of cell weights " << cWeights.size()
                 << " does not equal number of cells " << mesh_.nCells()
@@ -762,7 +763,7 @@ Foam::labelList Foam::parMetisDecomp::decompose
     Field<int> adjncy;
     // Offsets into adjncy
     Field<int> xadj;
-    metisDecomp::calcMetisCSR(globalCellCells, adjncy, xadj);
+    scotchDecomp::calcCSR(globalCellCells, adjncy, xadj);
 
     // decomposition options. 0 = use defaults
     List<int> options(3, 0);
diff --git a/src/surfMesh/Make/files b/src/surfMesh/Make/files
index 55b12af52e..a62f9c27f9 100644
--- a/src/surfMesh/Make/files
+++ b/src/surfMesh/Make/files
@@ -24,7 +24,6 @@ $(surfaceFormats)/ac3d/AC3DsurfaceFormatCore.C
 $(surfaceFormats)/ac3d/AC3DsurfaceFormatRunTime.C
 $(surfaceFormats)/ftr/FTRsurfaceFormatRunTime.C
 $(surfaceFormats)/gts/GTSsurfaceFormatRunTime.C
-$(surfaceFormats)/nas/NASsurfaceFormatCore.C
 $(surfaceFormats)/nas/NASsurfaceFormatRunTime.C
 $(surfaceFormats)/obj/OBJsurfaceFormatRunTime.C
 $(surfaceFormats)/off/OFFsurfaceFormatRunTime.C
diff --git a/src/surfMesh/Make/options b/src/surfMesh/Make/options
index e69de29bb2..7e207d0dbe 100644
--- a/src/surfMesh/Make/options
+++ b/src/surfMesh/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude
+
+LIB_LIBS = \
+    -lfileFormats
diff --git a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.H b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.H
index bf57a568a5..3708c213bc 100644
--- a/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.H
+++ b/src/surfMesh/surfaceFormats/nas/NASsurfaceFormat.H
@@ -47,7 +47,7 @@ SourceFiles
 #include "MeshedSurface.H"
 #include "MeshedSurfaceProxy.H"
 #include "UnsortedMeshedSurface.H"
-#include "NASsurfaceFormatCore.H"
+#include "NASCore.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,7 +64,7 @@ template<class Face>
 class NASsurfaceFormat
 :
     public MeshedSurface<Face>,
-    public NASsurfaceFormatCore
+    public NASCore
 {
     // Private Member Functions
 
diff --git a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.C b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.C
index 931088ea33..f85025ec4a 100644
--- a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.C
+++ b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.C
@@ -31,63 +31,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-bool Foam::fileFormats::STARCDsurfaceFormatCore::readHeader
-(
-    IFstream& is,
-    const word& signature
-)
-{
-    if (!is.good())
-    {
-        FatalErrorIn
-        (
-            "fileFormats::STARCDsurfaceFormatCore::readHeader(...)"
-        )
-            << "cannot read " << signature  << "  " << is.name()
-            << abort(FatalError);
-    }
-
-    word header;
-    label majorVersion;
-
-    string line;
-
-    is.getLine(line);
-    IStringStream(line)() >> header;
-
-    is.getLine(line);
-    IStringStream(line)() >> majorVersion;
-
-    // add other checks ...
-    if (header != signature)
-    {
-        Info<< "header mismatch " << signature << "  " << is.name()
-            << endl;
-    }
-
-    return true;
-}
-
-
-void Foam::fileFormats::STARCDsurfaceFormatCore::writeHeader
-(
-    Ostream& os,
-    const char* filetype
-)
-{
-    os  << "PROSTAR_" << filetype << nl
-        << 4000
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << " " << 0
-        << endl;
-}
-
-
 // parse things like this:
 //     CTNAME  1  someName
 // don't bother with the older comma-delimited format
@@ -102,12 +45,12 @@ Foam::fileFormats::STARCDsurfaceFormatCore::readInpCellTable
 
     regExp ctnameRE
     (
-        " *CTNA[^ ]*"       // keyword - min 4 chars
-        "[[:space:]]+"      // space delimited
-        "([0-9]+)"          // 1: <digits>
-        "[[:space:]]+"      // space delimited
+        " *CTNA[^ ]*"        // keyword - min 4 chars
+        "[[:space:]]+"       // space delimited
+        "([0-9]+)"           // 1: <digits>
+        "[[:space:]]+"       // space delimited
         "([^,[:space:]].*)", // 2: <name>
-        true                // ignore case
+        true                 // ignore case
     );
 
     string line;
@@ -132,78 +75,6 @@ Foam::fileFormats::STARCDsurfaceFormatCore::readInpCellTable
 }
 
 
-bool Foam::fileFormats::STARCDsurfaceFormatCore::readPoints
-(
-    IFstream& is,
-    pointField& points,
-    labelList& ids
-)
-{
-    //
-    // read .vrt file
-    // ~~~~~~~~~~~~~~
-
-    if (!is.good())
-    {
-        FatalErrorIn
-        (
-            "fileFormats::STARCDsurfaceFormatCore::readPoints(...)"
-        )
-            << "Cannot read file " << is.name()
-            << exit(FatalError);
-    }
-
-    readHeader(is, "PROSTAR_VERTEX");
-
-    DynamicList<point> dynPoints;
-    // STAR-CD index of points
-    DynamicList<label> dynPointId;
-
-    label lineLabel;
-    while ((is >> lineLabel).good())
-    {
-        scalar x, y, z;
-
-        is >> x >> y >> z;
-
-        dynPoints.append(point(x, y, z));
-        dynPointId.append(lineLabel);
-    }
-
-    points.transfer(dynPoints);
-    ids.transfer(dynPointId);
-
-    return true;
-}
-
-
-
-void Foam::fileFormats::STARCDsurfaceFormatCore::writePoints
-(
-    Ostream& os,
-    const pointField& pointLst
-)
-{
-    writeHeader(os, "VERTEX");
-
-    // Set the precision of the points data to 10
-    os.precision(10);
-
-    // force decimal point for Fortran input
-    os.setf(std::ios::showpoint);
-
-    forAll(pointLst, ptI)
-    {
-        os
-            << ptI + 1 << " "
-            << pointLst[ptI].x() << " "
-            << pointLst[ptI].y() << " "
-            << pointLst[ptI].z() << nl;
-    }
-    os.flush();
-}
-
-
 void Foam::fileFormats::STARCDsurfaceFormatCore::writeCase
 (
     Ostream& os,
@@ -238,4 +109,3 @@ void Foam::fileFormats::STARCDsurfaceFormatCore::writeCase
 
 
 // ************************************************************************* //
-
diff --git a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.H b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.H
index 0376d6248e..d99867d180 100644
--- a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.H
+++ b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormatCore.H
@@ -40,6 +40,7 @@ SourceFiles
 #include "Ostream.H"
 #include "OFstream.H"
 #include "MeshedSurface.H"
+#include "STARCDCore.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,21 +54,15 @@ namespace fileFormats
 \*---------------------------------------------------------------------------*/
 
 class STARCDsurfaceFormatCore
+:
+    public STARCDCore
 {
 protected:
 
     // Protected Member Functions
 
-    static bool readHeader(IFstream&, const word&);
-
-    static void writeHeader(Ostream&, const char* filetype);
-
     static Map<word> readInpCellTable(IFstream&);
 
-    static bool readPoints(IFstream&, pointField&, labelList& ids);
-
-    static void writePoints(Ostream&, const pointField&);
-
     static void writeCase
     (
         Ostream&,
diff --git a/src/triSurface/Make/options b/src/triSurface/Make/options
index 41306609f2..7e207d0dbe 100644
--- a/src/triSurface/Make/options
+++ b/src/triSurface/Make/options
@@ -1 +1,5 @@
-EXE_INC =
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude
+
+LIB_LIBS = \
+    -lfileFormats
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H
index 3438b60256..8c7014a900 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H
@@ -36,7 +36,7 @@ Description
     Example usage:
         myInterfacePatchName
         {
-            type                turbulentTemperatureCoupledBaffle;
+            type    compressible::turbulentTemperatureCoupledBaffle;
             neighbourFieldName  T;
             K                   K; // or none
             value               uniform 300;
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H
index 22df0cad0a..c4a6a7a763 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffleMixed/turbulentTemperatureCoupledBaffleMixedFvPatchScalarField.H
@@ -37,7 +37,7 @@ Description
     Example usage:
         myInterfacePatchName
         {
-            type                turbulentTemperatureCoupledBaffleMixed;
+            type        compressible::turbulentTemperatureCoupledBaffleMixed;
             neighbourFieldName  T;
             K                   K; // or none
             value               uniform 300;
diff --git a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
index 714a95ebbf..94282655bb 100644
--- a/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
+++ b/tutorials/heatTransfer/buoyantBoussinesqSimpleFoam/iglooWithFridges/system/snappyHexMeshDict
@@ -220,6 +220,12 @@ castellatedMeshControls
     // NOTE: This point should never be on a face, always inside a cell, even
     // after refinement.
     locationInMesh (3 0.28 0.43);
+
+
+    // Whether any faceZones (as specified in the refinementSurfaces)
+    // are only on the boundary of corresponding cellZones or also allow
+    // free-standing zone faces. Not used if there are no faceZones.
+    allowFreeStandingZoneFaces true;
 }
 
 
@@ -311,7 +317,8 @@ addLayersControls
     maxThicknessToMedialRatio 0.3;
 
     // Angle used to pick up medial axis points
-    minMedianAxisAngle 130;
+    // Note: changed(corrected) w.r.t 16x! 90 degrees corresponds to 130 in 16x.
+    minMedianAxisAngle 90;
 
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
index 70593293f5..948273e71d 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/snappyHexMeshDict
@@ -145,6 +145,12 @@ castellatedMeshControls
     // NOTE: This point should never be on a face, always inside a cell, even
     // after refinement.
     locationInMesh (3 3 0.43);
+
+
+    // Whether any faceZones (as specified in the refinementSurfaces)
+    // are only on the boundary of corresponding cellZones or also allow
+    // free-standing zone faces. Not used if there are no faceZones.
+    allowFreeStandingZoneFaces true;
 }
 
 
@@ -500,7 +506,9 @@ addLayersControls
     maxThicknessToMedialRatio 0.3;
 
     // Angle used to pick up medial axis points
-    minMedianAxisAngle 130;
+    // Note: changed(corrected) w.r.t 16x! 90 degrees corresponds to 130 in 16x.
+    minMedianAxisAngle 90;
+
 
     // Create buffer region for new layer terminations
     nBufferCellsNoExtrude 0;
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/U b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/U
index d692d79623..be47ccdfa5 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/U
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/U
@@ -33,17 +33,18 @@ boundaryField
     {
         type            flowRateInletVelocity;
         flowRate        0.00379;
-        value           uniform (-0 10.82857143 -0);
+        value           uniform (0 14.68 0);
     }
     inletSides
     {
         type            flowRateInletVelocity;
         flowRate        0.00832;
-        value           uniform (-0 11.88571429 -0);
+        value           uniform (0 17.79 0);
     }
     outlet
     {
-        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
     }
     walls
     {
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/k b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/k
index 59479eef3a..d31f55adb4 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/k
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/k
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 2 -2 0 0 0 0];
 
-internalField   uniform 3.75e-11;
+internalField   uniform 3.75e-9;
 
 boundaryField
 {
@@ -33,22 +33,23 @@ boundaryField
     {
         type            turbulentIntensityKineticEnergyInlet;
         intensity       0.15;
-        value           uniform 3.75e-11;
+        value           uniform 3.75e-9;
     }
     inletSides
     {
         type            turbulentIntensityKineticEnergyInlet;
         intensity       0.16;
-        value           uniform 3.75e-11;
+        value           uniform 3.75e-9;
     }
     outlet
     {
-        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform 3.75e-9;
     }
     walls
     {
         type            compressible::kqRWallFunction;
-        value           uniform 3.75e-11;
+        value           uniform 0;
     }
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/omega b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/omega
index f7c68ea1c4..0538e96fb7 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/omega
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0.org/omega
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 0 -1 0 0 0 0];
 
-internalField   uniform 0.0016;
+internalField   uniform 4.5e-3;
 
 boundaryField
 {
@@ -34,18 +34,19 @@ boundaryField
         type            compressible::turbulentMixingLengthFrequencyInlet;
         mixingLength    0.007;
         k               k;
-        value           uniform 0.0016;
+        value           uniform 4.5e-3;
     }
     inletSides
     {
         type            compressible::turbulentMixingLengthFrequencyInlet;
         mixingLength    0.007;
         k               k;
-        value           uniform 0.0016;
+        value           uniform 4.5e-3;
     }
     outlet
     {
-        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform 4.5e-3;
     }
     walls
     {
@@ -53,7 +54,7 @@ boundaryField
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.0016;
+        value           uniform 0;
     }
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/U b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/U
index d692d79623..b2ed1ea6d7 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/U
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/U
@@ -43,7 +43,8 @@ boundaryField
     }
     outlet
     {
-        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform (0 0 0);
     }
     walls
     {
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/k b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/k
index 59479eef3a..d31f55adb4 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/k
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/k
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 2 -2 0 0 0 0];
 
-internalField   uniform 3.75e-11;
+internalField   uniform 3.75e-9;
 
 boundaryField
 {
@@ -33,22 +33,23 @@ boundaryField
     {
         type            turbulentIntensityKineticEnergyInlet;
         intensity       0.15;
-        value           uniform 3.75e-11;
+        value           uniform 3.75e-9;
     }
     inletSides
     {
         type            turbulentIntensityKineticEnergyInlet;
         intensity       0.16;
-        value           uniform 3.75e-11;
+        value           uniform 3.75e-9;
     }
     outlet
     {
-        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform 3.75e-9;
     }
     walls
     {
         type            compressible::kqRWallFunction;
-        value           uniform 3.75e-11;
+        value           uniform 0;
     }
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/omega b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/omega
index f7c68ea1c4..f9b551e72b 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/omega
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/0/omega
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 0 -1 0 0 0 0];
 
-internalField   uniform 0.0016;
+internalField   uniform 4.5e-3;
 
 boundaryField
 {
@@ -34,18 +34,19 @@ boundaryField
         type            compressible::turbulentMixingLengthFrequencyInlet;
         mixingLength    0.007;
         k               k;
-        value           uniform 0.0016;
+        value           uniform 4.5e-3;
     }
     inletSides
     {
         type            compressible::turbulentMixingLengthFrequencyInlet;
         mixingLength    0.007;
         k               k;
-        value           uniform 0.0016;
+        value           uniform 4.5-3;
     }
     outlet
     {
-        type            zeroGradient;
+        type            inletOutlet;
+        inletValue      uniform 4.5e-3;
     }
     walls
     {
@@ -53,7 +54,7 @@ boundaryField
         Cmu             0.09;
         kappa           0.41;
         E               9.8;
-        value           uniform 0.0016;
+        value           uniform 0;
     }
 }
 
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allrun b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allrun
index 7aca4e2d34..6eb3715e65 100755
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allrun
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/Allrun
@@ -5,5 +5,9 @@
 # create mesh
 runApplication blockMesh
 
+# initialise with potentialFoam solution
+runApplication potentialFoam
+rm -f 0/phi
+
 # run the solver
 runApplication porousExplicitSourceReactingParcelFoam
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/g b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/g
index 51944e7abd..27d4d32488 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/g
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/g
@@ -16,7 +16,7 @@ FoamFile
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 dimensions      [0 1 -2 0 0 0 0];
-value           ( 0 0 0 );
+value           ( 0 -9.81 0 );
 
 
 // ************************************************************************* //
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
index f62411e787..e8088bf7bd 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/polyMesh/boundary
@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
diff --git a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
index d9e40ac14d..ddb7eaf3a9 100644
--- a/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
+++ b/tutorials/lagrangian/porousExplicitSourceReactingParcelFoam/verticalChannel/constant/reactingCloud1Properties
@@ -129,7 +129,7 @@ StandardWallInteractionCoeffs
 
 RanzMarshallCoeffs
 {
-    BirdCorrection  on; // off;
+    BirdCorrection  off;
 }
 
 SinglePhaseMixtureCoeffs