diff --git a/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H b/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H
deleted file mode 100644
index d024bd2017..0000000000
--- a/applications/solvers/discreteMethods/dsmc/dsmcFoam/createFields.H
+++ /dev/null
@@ -1,162 +0,0 @@
-
-    Info<< nl << "Reading field boundaryT" << endl;
-    volScalarField boundaryT
-    (
-        IOobject
-        (
-            "boundaryT",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field boundaryU" << endl;
-    volVectorField boundaryU
-    (
-        IOobject
-        (
-            "boundaryU",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field rhoN (number density)" << endl;
-    volScalarField rhoN
-    (
-        IOobject
-        (
-            "rhoN",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field rhoM (mass density)" << endl;
-    volScalarField rhoM
-    (
-        IOobject
-        (
-            "rhoM",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field rhoNdsmc (dsmc particle density)" << endl;
-    volScalarField dsmcRhoN
-    (
-        IOobject
-        (
-            "dsmcRhoN",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field momentum (momentum density)" << endl;
-    volVectorField momentum
-    (
-        IOobject
-        (
-            "momentum",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field linearKE (linear kinetic energy density)"
-        << endl;
-
-    volScalarField linearKE
-    (
-        IOobject
-        (
-            "linearKE",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field internalE (internal energy density)" << endl;
-    volScalarField internalE
-    (
-        IOobject
-        (
-            "internalE",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field iDof (internal degree of freedom density)"
-        << endl;
-
-    volScalarField iDof
-    (
-        IOobject
-        (
-            "iDof",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field q (surface heat transfer)" << endl;
-    volScalarField q
-    (
-        IOobject
-        (
-            "q",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Reading field fD (surface force density)" << endl;
-    volVectorField fD
-    (
-        IOobject
-        (
-            "fD",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< nl << "Constructing dsmcCloud " << endl;
-
-    dsmcCloud dsmc("dsmc", boundaryT, boundaryU);
diff --git a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
index 84b499451f..f5a57f15e8 100644
--- a/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
+++ b/applications/solvers/discreteMethods/dsmc/dsmcFoam/dsmcFoam.C
@@ -41,53 +41,21 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "createFields.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+    Info<< nl << "Constructing dsmcCloud " << endl;
+
+    dsmcCloud dsmc("dsmc", mesh);
+
     Info<< "\nStarting time loop\n" << endl;
 
-    while (runTime.run())
+    while (runTime.loop())
     {
-        runTime++;
-
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        // Carry out dsmcCloud timestep
-
         dsmc.evolve();
 
-        // Retrieve flow field data from dsmcCloud
-
-        rhoN = dsmc.rhoN();
-        rhoN.correctBoundaryConditions();
-
-        rhoM = dsmc.rhoM();
-        rhoM.correctBoundaryConditions();
-
-        dsmcRhoN = dsmc.dsmcRhoN();
-        dsmcRhoN.correctBoundaryConditions();
-
-        momentum = dsmc.momentum();
-        momentum.correctBoundaryConditions();
-
-        linearKE = dsmc.linearKE();
-        linearKE.correctBoundaryConditions();
-
-        internalE = dsmc.internalE();
-        internalE.correctBoundaryConditions();
-
-        iDof = dsmc.iDof();
-        iDof.correctBoundaryConditions();
-
-        // Retrieve surface field data from dsmcCloud
-
-        q = dsmc.q();
-
-        fD = dsmc.fD();
-
-        // Print status of dsmcCloud
-
         dsmc.info();
 
         runTime.write();
diff --git a/applications/test/DimensionedField/DimensionedFieldTest.C b/applications/test/DimensionedField/DimensionedFieldTest.C
deleted file mode 100644
index bc1f875dbd..0000000000
--- a/applications/test/DimensionedField/DimensionedFieldTest.C
+++ /dev/null
@@ -1,133 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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
-
-\*---------------------------------------------------------------------------*/
-
-#include "argList.H"
-#include "Time.H"
-#include "DimensionedFields.H"
-#include "DimensionedSphericalTensorField.H"
-#include "vector.H"
-#include "tensor.H"
-#include "GeoMesh.H"
-
-using namespace Foam;
-
-namespace Foam
-{
-
-class vMesh
-{
-
-public:
-
-    vMesh()
-    {}
-
-    label size() const
-    {
-        return 10;
-    }
-};
-
-};
-
-template<>
-const word Foam::DimensionedField<scalar, GeoMesh<vMesh> >::typeName
-(
-    "dimenionedScalarField"
-);
-
-template<>
-const word Foam::DimensionedField<vector, GeoMesh<vMesh> >::typeName
-(
-    "dimenionedVectorField"
-);
-
-template<>
-const word Foam::DimensionedField<tensor, GeoMesh<vMesh> >::typeName
-(
-    "dimenionedTensorField"
-);
-
-template<>
-const word Foam::DimensionedField<sphericalTensor, GeoMesh<vMesh> >::typeName
-(
-    "dimenionedSphericalTensorField"
-);
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Main program:
-
-int main(int argc, char *argv[])
-{
-
-#   include "setRootCase.H"
-#   include "createTime.H"
-
-    vMesh vm;
-
-    DimensionedField<scalar, GeoMesh<vMesh> > dsf
-    (
-        IOobject
-        (
-            "dsf",
-            runTime.timeName(),
-            runTime,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        vm
-    );
-
-    Info<< dsf << endl;
-    dsf += dsf;
-    dsf -= dimensionedScalar("5", dsf.dimensions(), 5.0);
-    Info<< dsf << endl;
-
-    Info<< sqr(dsf + dsf) - sqr(dsf + dsf) << endl;
-
-    DimensionedField<vector, GeoMesh<vMesh> > dvf
-    (
-        IOobject
-        (
-            "dvf",
-            runTime.timeName(),
-            runTime,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        vm
-    );
-
-    Info<< (dvf ^ (dvf ^ dvf)) << endl;
-
-    Info << "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //
diff --git a/applications/test/DimensionedField/Make/files b/applications/test/DimensionedField/Make/files
deleted file mode 100644
index 02c479015e..0000000000
--- a/applications/test/DimensionedField/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-DimensionedFieldTest.C
-
-EXE = $(FOAM_USER_APPBIN)/DimensionedFieldTest
diff --git a/applications/test/DimensionedField/Make/options b/applications/test/DimensionedField/Make/options
deleted file mode 100644
index 6a9e9810b3..0000000000
--- a/applications/test/DimensionedField/Make/options
+++ /dev/null
@@ -1,2 +0,0 @@
-/* EXE_INC = -I$(LIB_SRC)/cfdTools/include */
-/* EXE_LIBS = -lfiniteVolume */
diff --git a/applications/test/HashTable3/hashTableTest3.C b/applications/test/HashTable3/hashTableTest3.C
index 910e4f2290..6229e6a304 100644
--- a/applications/test/HashTable3/hashTableTest3.C
+++ b/applications/test/HashTable3/hashTableTest3.C
@@ -33,7 +33,6 @@ Description
 #include "HashPtrTable.H"
 #include "Map.H"
 #include "StaticHashTable.H"
-#include "HashTbl.H"
 #include "cpuTime.H"
 
 using namespace Foam;
@@ -53,7 +52,7 @@ int main(int argc, char *argv[])
     // Map<label> map(2 * nSize);
     // HashTable<label, label, Hash<label> > map(2 * nSize);
     // StaticHashTable<label, label, Hash<label> > map(2 * nSize);
-    HashTbl<label, label, Hash<label> > map(2 * nSize);
+    HashTable<label, label, Hash<label> > map(2 * nSize);
 
     Info<< "Constructed map of size: " << nSize
         << " (size " << map.size() << " capacity " << map.capacity() << ") "
diff --git a/applications/test/Matrix/MatrixTest.C b/applications/test/Matrix/MatrixTest.C
index 24bcc68257..5576291fd5 100644
--- a/applications/test/Matrix/MatrixTest.C
+++ b/applications/test/Matrix/MatrixTest.C
@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Matrix.H"
+#include "SquareMatrix.H"
 #include "vector.H"
 
 using namespace Foam;
@@ -34,7 +34,7 @@ using namespace Foam;
 
 int main(int argc, char *argv[])
 {
-    Matrix<scalar> hmm(3, 3);
+    SquareMatrix<scalar> hmm(3);
 
     hmm[0][0] = -3.0;
     hmm[0][1] = 10.0;
@@ -46,27 +46,27 @@ int main(int argc, char *argv[])
     hmm[2][1] = 6.0;
     hmm[2][2] = 1.0;
 
-    Info<< hmm << endl << hmm - 2.0*(-hmm) << endl;
+    //Info<< hmm << endl << hmm - 2.0*(-hmm) << endl;
     Info<< max(hmm) << endl;
     Info<< min(hmm) << endl;
 
-    Matrix<scalar> hmm2(3, 3, 1.0);
+    SquareMatrix<scalar> hmm2(3, 1.0);
 
     hmm = hmm2;
 
     Info<< hmm << endl;
 
-    Matrix<scalar> hmm3(Sin);
+    SquareMatrix<scalar> hmm3(Sin);
 
     Info<< hmm3 << endl;
 
-    Matrix<scalar> hmm4;
+    SquareMatrix<scalar> hmm4;
 
     hmm4 = hmm2;
 
     Info<< hmm4 << endl;
 
-    Matrix<scalar> hmm5;
+    SquareMatrix<scalar> hmm5;
 
     hmm4 = hmm5;
     Info<< hmm5 << endl;
diff --git a/applications/test/ODETest/ODETest.C b/applications/test/ODETest/ODETest.C
index dfc7a663bb..2d8cecfa9a 100644
--- a/applications/test/ODETest/ODETest.C
+++ b/applications/test/ODETest/ODETest.C
@@ -69,7 +69,7 @@ public:
         const scalar x,
         const scalarField& y,
         scalarField& dfdx,
-        Matrix<scalar>& dfdy
+        scalarSquareMatrix& dfdy
     ) const
     {
         dfdx[0] = 0.0;
diff --git a/applications/test/PtrList/PtrListTest.C b/applications/test/PtrList/PtrListTest.C
index 5413a662b8..7fbff08d11 100644
--- a/applications/test/PtrList/PtrListTest.C
+++ b/applications/test/PtrList/PtrListTest.C
@@ -57,6 +57,8 @@ public:
         Info <<"delete Scalar: " << data_ << endl;
     }
 
+    autoPtr<Scalar> clone() const;
+
     friend Ostream& operator<<(Ostream& os, const Scalar& val)
     {
         os << val.data_;
diff --git a/applications/test/UIndirectListTest/UIndirectListTest.C b/applications/test/UIndirectListTest/UIndirectListTest.C
index 6902f47eea..67943b7d94 100644
--- a/applications/test/UIndirectListTest/UIndirectListTest.C
+++ b/applications/test/UIndirectListTest/UIndirectListTest.C
@@ -84,11 +84,11 @@ int main(int argc, char *argv[])
 
     // test List operations
 
-    List<double> flatList = UIndirectList<double>(completeList, addresses);
-    Info<< "List assigned from UIndirectList: " << flatList << endl;
+    List<double> flatList(UIndirectList<double>(completeList, addresses));
+    Info<< "List constructed from UIndirectList: " << flatList << endl;
 
-    List<double> flatList2(UIndirectList<double>(completeList, addresses));
-    Info<< "List constructed from UIndirectList: " << flatList2 << endl;
+    flatList = UIndirectList<double>(completeList, addresses);
+    Info<< "List assigned from UIndirectList: " << flatList << endl;
 
     flatList.append(UIndirectList<double>(completeList, addresses));
     Info<< "List::append(UIndirectList): " << flatList << endl;
diff --git a/applications/test/parallel/parallelTest.C b/applications/test/parallel/parallelTest.C
index 6120de31a5..caf9e6fc8d 100644
--- a/applications/test/parallel/parallelTest.C
+++ b/applications/test/parallel/parallelTest.C
@@ -84,7 +84,7 @@ int main(int argc, char *argv[])
         // Sync how many to send
         labelListList allNTrans(Pstream::nProcs());
         allNTrans[Pstream::myProcNo()] = nSend;
-        combineReduce(allNTrans, mapDistribute::listEq());
+        combineReduce(allNTrans, UPstream::listEq());
 
         // Collect items to be sent
         labelListList sendMap(Pstream::nProcs());
@@ -161,7 +161,7 @@ int main(int argc, char *argv[])
                 toMaster << data;
             }
 
-            Perr<< "slave receiving from master " 
+            Perr<< "slave receiving from master "
                 << Pstream::masterNo() << endl;
             IPstream fromMaster(Pstream::blocking, Pstream::masterNo());
             fromMaster >> data;
diff --git a/applications/test/slicedField/slicedFieldTest.C b/applications/test/slicedField/slicedFieldTest.C
index a82b4ac712..dd8da25223 100644
--- a/applications/test/slicedField/slicedFieldTest.C
+++ b/applications/test/slicedField/slicedFieldTest.C
@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    volField
+    slicedFieldTest
 
 Description
 
@@ -32,16 +32,16 @@ Description
 #include "fvCFD.H"
 #include "SlicedGeometricField.H"
 #include "slicedFvPatchFields.H"
+#include "slicedSurfaceFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.H"
+    #include "createTime.H"
+    #include "createMesh.H"
 
     Info<< "Reading field p\n" << endl;
     volScalarField p
@@ -91,7 +91,7 @@ int main(int argc, char *argv[])
     Info<< C << endl;
     Info<< (C & U) << endl;
 
-    SlicedGeometricField<vector, fvPatchField, slicedFvPatchField, surfaceMesh>
+    SlicedGeometricField<vector, fvsPatchField, slicedFvsPatchField, surfaceMesh>
     Sf
     (
         IOobject
@@ -105,7 +105,7 @@ int main(int argc, char *argv[])
         mesh.faceAreas()
     );
 
-    Info<< Sf << endl;
+    //Info<< Sf << endl;
 
     return 0;
 }
diff --git a/applications/test/testPointEdgeWave/testPointEdgeWave.C b/applications/test/testPointEdgeWave/testPointEdgeWave.C
index 3a1ff900eb..124b7a921a 100644
--- a/applications/test/testPointEdgeWave/testPointEdgeWave.C
+++ b/applications/test/testPointEdgeWave/testPointEdgeWave.C
@@ -48,8 +48,6 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createPolyMesh.H"
 
-    pointMesh pMesh(mesh);
-
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     // Get name of patch
@@ -99,7 +97,7 @@ int main(int argc, char *argv[])
 
     PointEdgeWave<pointEdgePoint> wallCalc
     (
-        pMesh,
+        mesh,
         wallPoints,
         wallInfo,
 
@@ -119,7 +117,7 @@ int main(int argc, char *argv[])
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        pMesh,
+        pointMesh::New(mesh),
         dimensionedScalar("wallDist", dimLength, 0.0)
     );
 
diff --git a/applications/test/volField/volFieldTest.C b/applications/test/volField/volFieldTest.C
index 3b7dafef8e..f910c22544 100644
--- a/applications/test/volField/volFieldTest.C
+++ b/applications/test/volField/volFieldTest.C
@@ -34,11 +34,10 @@ Application
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
 
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.H"
+    #include "createTime.H"
+    #include "createMesh.H"
 
     Info<< "Reading field p\n" << endl;
     volScalarField p
@@ -70,12 +69,9 @@ int main(int argc, char *argv[])
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
-    //Info<< transform(dimensionedTensor("I", dimless, 0.1*I), U) << endl;
-
-
-    GeometricField<sphericalTensor, fvPatchField, volMesh> st
+    GeometricField<symmTensor, fvPatchField, volMesh> st
     (
         IOobject
         (
@@ -86,8 +82,8 @@ int main(int argc, char *argv[])
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensioned<sphericalTensor>("st", dimless, sphericalTensor::I),
-        zeroGradientFvPatchSphericalTensorField::typeName
+        dimensioned<symmTensor>("st", dimless, symmTensor::one),
+        zeroGradientFvPatchSymmTensorField::typeName
     );
 
     //Info<< fvc::div(st) << endl;
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index be267f8724..b61d997399 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -45,7 +45,7 @@ Description
     so cannot have any cells in any other zone.
     - useCellZonesOnly does not do a walk and uses the cellZones only. Use
     this if you don't mind having disconnected domains in a single region.
-    This option requires all cells to be in one (and one only) region.
+    This option requires all cells to be in one (and one only) cellZone.
 
 \*---------------------------------------------------------------------------*/
 
diff --git a/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C
index 95077df225..34f231850f 100644
--- a/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C
+++ b/applications/utilities/postProcessing/miscellaneous/dsmcFieldsCalc/dsmcFieldsCalc.C
@@ -113,7 +113,16 @@ void Foam::calc(const argList& args, const Time& runTime, const fvMesh& mesh)
         return;
     }
 
-    wordList extensiveVVFNames(IStringStream ("(momentumMean)")());
+    wordList extensiveVVFNames
+    (
+        IStringStream
+        (
+            "( \
+                momentumMean \
+                fDMean \
+            )"
+        )()
+    );
 
     PtrList<volVectorField> extensiveVVFs(extensiveVVFNames.size());
 
diff --git a/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
index f8b831001d..9e7725ec23 100644
--- a/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
+++ b/applications/utilities/preProcessing/dsmcInitialise/dsmcInitialise.C
@@ -44,9 +44,21 @@ int main(int argc, char *argv[])
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+    IOdictionary dsmcInitialiseDict
+    (
+        IOobject
+        (
+            "dsmcInitialiseDict",
+            mesh.time().system(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
     Info<< "Initialising dsmc for Time = " << runTime.timeName() << nl << endl;
 
-    dsmcCloud dsmc("dsmc", mesh);
+    dsmcCloud dsmc("dsmc", mesh, dsmcInitialiseDict);
 
     label totalMolecules = dsmc.size();
 
diff --git a/bin/mpirunDebug b/bin/mpirunDebug
index 57df760873..2376d22601 100755
--- a/bin/mpirunDebug
+++ b/bin/mpirunDebug
@@ -206,6 +206,21 @@ done
 cmd=""
 if [ .$WM_MPLIB = .OPENMPI ]; then
     cmd="mpirun -app $PWD/mpirun.schema </dev/null"
+elif [ .$WM_MPLIB = .MPICH ]; then
+    cmd="mpiexec"
+    for ((proc=0; proc<$nProcs; proc++))
+    do
+        read procCmd
+
+        procXtermCmdFile="$PWD/processor${proc}Xterm.sh"
+        echo "#!/bin/sh" > $procXtermCmdFile
+        echo "$procCmd" >> $procXtermCmdFile
+        chmod +x $procXtermCmdFile
+        if [ $proc -ne 0 ]; then
+            cmd="${cmd} :"
+        fi
+        cmd="${cmd} -n 1 ${procXtermCmdFile}"
+    done < $PWD/mpirun.schema
 fi
 
 echo "Constructed $PWD/mpirun.schema file."
diff --git a/etc/controlDict b/etc/controlDict
index 97af349e7a..ef72407eaf 100644
--- a/etc/controlDict
+++ b/etc/controlDict
@@ -377,6 +377,7 @@ DebugSwitches
     displacementLaplacian 0;
     displacementSBRStress 0;
     distanceSurface     0;
+    distribution        0;
     downwind            0;
     dragModel           0;
     duplicatePoints     0;
diff --git a/src/OpenFOAM/containers/HashTables/HashSet/HashSet.C b/src/OpenFOAM/containers/HashTables/HashSet/HashSet.C
index d63921f95b..4701c8020d 100644
--- a/src/OpenFOAM/containers/HashTables/HashSet/HashSet.C
+++ b/src/OpenFOAM/containers/HashTables/HashSet/HashSet.C
@@ -105,7 +105,7 @@ template<class Key, class Hash>
 void Foam::HashSet<Key, Hash>::operator&=(const HashSet<Key, Hash>& rhs)
 {
     // Remove elements not also found in rhs
-    for (iterator iter = this->cbegin(); iter != this->cend(); ++iter)
+    for (iterator iter = this->begin(); iter != this->end(); ++iter)
     {
         if (!rhs.found(iter.key()))
         {
@@ -145,8 +145,6 @@ void Foam::HashSet<Key, Hash>::operator-=(const HashSet<Key, Hash>& rhs)
 }
 
 
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
 /* * * * * * * * * * * * * * * * Global operators  * * * * * * * * * * * * * */
 
 template<class Key, class Hash>
diff --git a/src/OpenFOAM/containers/Lists/UList/UList.H b/src/OpenFOAM/containers/Lists/UList/UList.H
index 087f685ff6..b829c99ea5 100644
--- a/src/OpenFOAM/containers/Lists/UList/UList.H
+++ b/src/OpenFOAM/containers/Lists/UList/UList.H
@@ -184,6 +184,12 @@ public:
         //  an out-of-range element returns false without any ill-effects
         inline const T& operator[](const label) const;
 
+        //- Return last element of UList.
+        inline T& last();
+
+        //- Return last element of UList.
+        inline const T& last() const;
+
         //- Allow cast to a const List<T>&
         inline operator const Foam::List<T>&() const;
 
diff --git a/src/OpenFOAM/containers/Lists/UList/UListI.H b/src/OpenFOAM/containers/Lists/UList/UListI.H
index 5c1df4be15..aa7d396f10 100644
--- a/src/OpenFOAM/containers/Lists/UList/UListI.H
+++ b/src/OpenFOAM/containers/Lists/UList/UListI.H
@@ -114,6 +114,20 @@ inline void Foam::UList<T>::checkIndex(const label i) const
 }
 
 
+template<class T>
+inline T& Foam::UList<T>::last()
+{
+    return this->operator[](this->size()-1);
+}
+
+
+template<class T>
+inline const T& Foam::UList<T>::last() const
+{
+    return this->operator[](this->size()-1);
+}
+
+
 template<class T>
 inline const T* Foam::UList<T>::cdata() const
 {
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C
index 87f940db40..41cb1c4f11 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C
@@ -82,9 +82,6 @@ void Foam::PstreamBuffers::finishedSends(labelListList& sizes, const bool block)
 
     if (commsType_ == UPstream::nonBlocking)
     {
-        labelListList sizes;
-        labelListList send,recv;
-
         Pstream::exchange<DynamicList<char>, char>
         (
             sendBuf_,
@@ -96,20 +93,30 @@ void Foam::PstreamBuffers::finishedSends(labelListList& sizes, const bool block)
     }
     else
     {
-        sizes.setSize(UPstream::nProcs());
-        labelList& nsTransPs = sizes[UPstream::myProcNo()];
-        nsTransPs.setSize(UPstream::nProcs());
+        FatalErrorIn
+        (
+            "PstreamBuffers::finishedSends(labelListList&, const bool)"
+        )   << "Obtaining sizes not supported in "
+            << UPstream::commsTypeNames[commsType_] << endl
+            << " since transfers already in progress. Use non-blocking instead."
+            << exit(FatalError);
 
-        forAll(sendBuf_, procI)
-        {
-            nsTransPs[procI] = sendBuf_[procI].size();
-        }
-
-        // Send sizes across.
-        int oldTag = UPstream::msgType();
-        UPstream::msgType() = tag_;
-        combineReduce(sizes, UPstream::listEq());
-        UPstream::msgType() = oldTag;
+        // Note: possible only if using different tag from write started
+        // by ~UOPstream. Needs some work.
+        //sizes.setSize(UPstream::nProcs());
+        //labelList& nsTransPs = sizes[UPstream::myProcNo()];
+        //nsTransPs.setSize(UPstream::nProcs());
+        //
+        //forAll(sendBuf_, procI)
+        //{
+        //    nsTransPs[procI] = sendBuf_[procI].size();
+        //}
+        //
+        //// Send sizes across.
+        //int oldTag = UPstream::msgType();
+        //UPstream::msgType() = tag_;
+        //combineReduce(sizes, UPstream::listEq());
+        //UPstream::msgType() = oldTag;
     }
 }
 
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H
index d3dbbf23db..75f26f1d6a 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H
@@ -138,7 +138,8 @@ public:
         void finishedSends(const bool block = true);
 
         //- Mark all sends as having been done. Same as above but also returns
-        //  sizes (bytes) transferred.
+        //  sizes (bytes) transferred. Note:currently only valid for
+        //  non-blocking.
         void finishedSends(labelListList& sizes, const bool block = true);
 
 };
diff --git a/src/OpenFOAM/db/IOstreams/Pstreams/exchange.C b/src/OpenFOAM/db/IOstreams/Pstreams/exchange.C
index 64ced34ecf..b950ed22f0 100644
--- a/src/OpenFOAM/db/IOstreams/Pstreams/exchange.C
+++ b/src/OpenFOAM/db/IOstreams/Pstreams/exchange.C
@@ -50,43 +50,42 @@ void Pstream::exchange
     const bool block
 )
 {
-    if (UPstream::parRun())
+    if (!contiguous<T>())
     {
-        if (!contiguous<T>())
-        {
-            FatalErrorIn
-            (
-                "Pstream::exchange(..)"
-            )   << "Continuous data only." << Foam::abort(FatalError);
-        }
+        FatalErrorIn
+        (
+            "Pstream::exchange(..)"
+        )   << "Continuous data only." << Foam::abort(FatalError);
+    }
 
-        if (sendBufs.size() != UPstream::nProcs())
-        {
-            FatalErrorIn
-            (
-                "Pstream::exchange(..)"
-            )   << "Size of list:" << sendBufs.size()
-                << " does not equal the number of processors:"
-                << UPstream::nProcs()
-                << Foam::abort(FatalError);
-        }
+    if (sendBufs.size() != UPstream::nProcs())
+    {
+        FatalErrorIn
+        (
+            "Pstream::exchange(..)"
+        )   << "Size of list:" << sendBufs.size()
+            << " does not equal the number of processors:"
+            << UPstream::nProcs()
+            << Foam::abort(FatalError);
+    }
 
-        sizes.setSize(UPstream::nProcs());
-        labelList& nsTransPs = sizes[UPstream::myProcNo()];
-        nsTransPs.setSize(UPstream::nProcs());
+    sizes.setSize(UPstream::nProcs());
+    labelList& nsTransPs = sizes[UPstream::myProcNo()];
+    nsTransPs.setSize(UPstream::nProcs());
 
-        forAll(sendBufs, procI)
-        {
-            nsTransPs[procI] = sendBufs[procI].size();
-        }
-
-        // Send sizes across.
-        int oldTag = UPstream::msgType();
-        UPstream::msgType() = tag;
-        combineReduce(sizes, UPstream::listEq());
-        UPstream::msgType() = oldTag;
+    forAll(sendBufs, procI)
+    {
+        nsTransPs[procI] = sendBufs[procI].size();
+    }
 
+    // Send sizes across.
+    int oldTag = UPstream::msgType();
+    UPstream::msgType() = tag;
+    combineReduce(sizes, UPstream::listEq());
+    UPstream::msgType() = oldTag;
 
+    if (Pstream::parRun())
+    {
         // Set up receives
         // ~~~~~~~~~~~~~~~
 
diff --git a/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C b/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C
index 7fe2291357..f7040df07c 100644
--- a/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C
+++ b/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.C
@@ -37,8 +37,43 @@ void Foam::OutputFilterFunctionObject<OutputFilter>::readDict()
     dict_.readIfPresent("region", regionName_);
     dict_.readIfPresent("dictionary", dictName_);
     dict_.readIfPresent("enabled", enabled_);
+    dict_.readIfPresent("storeFilter", storeFilter_);
 }
 
+template<class OutputFilter>
+void Foam::OutputFilterFunctionObject<OutputFilter>::allocateFilter()
+{
+    if (dictName_.size())
+    {
+        ptr_.reset
+        (
+            new IOOutputFilter<OutputFilter>
+            (
+                name(),
+                time_.lookupObject<objectRegistry>(regionName_),
+                dictName_
+            )
+        );
+    }
+    else
+    {
+        ptr_.reset
+        (
+            new OutputFilter
+            (
+                name(),
+                time_.lookupObject<objectRegistry>(regionName_),
+                dict_
+            )
+        );
+    }
+}
+
+template<class OutputFilter>
+void Foam::OutputFilterFunctionObject<OutputFilter>::destroyFilter()
+{
+    ptr_.reset();
+}
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -56,6 +91,7 @@ Foam::OutputFilterFunctionObject<OutputFilter>::OutputFilterFunctionObject
     regionName_(polyMesh::defaultRegion),
     dictName_(),
     enabled_(true),
+    storeFilter_(true),
     outputControl_(t, dict)
 {
     readDict();
@@ -83,32 +119,9 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::start()
 {
     readDict();
 
-    if (enabled_)
+    if (enabled_&&storeFilter_)
     {
-        if (dictName_.size())
-        {
-            ptr_.reset
-            (
-                new IOOutputFilter<OutputFilter>
-                (
-                    name(),
-                    time_.lookupObject<objectRegistry>(regionName_),
-                    dictName_
-                )
-            );
-        }
-        else
-        {
-            ptr_.reset
-            (
-                new OutputFilter
-                (
-                    name(),
-                    time_.lookupObject<objectRegistry>(regionName_),
-                    dict_
-                )
-            );
-        }
+        allocateFilter();
     }
 
     return true;
@@ -120,12 +133,22 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::execute()
 {
     if (enabled_)
     {
+        if (!storeFilter_)
+        {
+            allocateFilter();
+        }
+
         ptr_->execute();
 
         if (outputControl_.output())
         {
             ptr_->write();
         }
+
+        if (!storeFilter_)
+        {
+            destroyFilter();
+        }
     }
 
     return true;
@@ -137,12 +160,22 @@ bool Foam::OutputFilterFunctionObject<OutputFilter>::end()
 {
     if (enabled_)
     {
+        if (!storeFilter_)
+        {
+            allocateFilter();
+        }
+
         ptr_->end();
 
         if (outputControl_.output())
         {
             ptr_->write();
         }
+
+        if (!storeFilter_)
+        {
+            destroyFilter();
+        }
     }
 
     return true;
diff --git a/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.H b/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.H
index f6c183e803..47009d2643 100644
--- a/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.H
+++ b/src/OpenFOAM/db/functionObjects/OutputFilterFunctionObject/OutputFilterFunctionObject.H
@@ -78,6 +78,10 @@ class OutputFilterFunctionObject
         //- Switch for the execution of the functionObject
         bool enabled_;
 
+        //- Switch to store filter in between writes or use on-the-fly
+        //  construction
+        bool storeFilter_;
+
         //- Output controls
         outputFilterOutputControl outputControl_;
 
@@ -89,6 +93,12 @@ class OutputFilterFunctionObject
 
         //- Read relevant dictionary entries
         void readDict();
+ 
+        //- Creates most of the data associated with this object.
+        void allocateFilter();
+
+        //- Destroys most of the data associated with this object.
+        void destroyFilter();
 
         //- Disallow default bitwise copy construct
         OutputFilterFunctionObject(const OutputFilterFunctionObject&);
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.C b/src/OpenFOAM/matrices/Matrix/Matrix.C
index 7f91a5f36c..0062d897a4 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.C
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.C
@@ -235,12 +235,12 @@ void Foam::Matrix<Form, Type>::operator=(const Matrix<Form, Type>& a)
 template<class Form, class Type>
 const Type& Foam::max(const Matrix<Form, Type>& a)
 {
-    label nm = a.n_*a.m_;
+    label nm = a.n()*a.m();
 
     if (nm)
     {
         label curMaxI = 0;
-        const Type* v = a.v_[0];
+        const Type* v = a[0];
 
         for (register label i=1; i<nm; i++)
         {
@@ -267,12 +267,12 @@ const Type& Foam::max(const Matrix<Form, Type>& a)
 template<class Form, class Type>
 const Type& Foam::min(const Matrix<Form, Type>& a)
 {
-    label nm = a.n_*a.m_;
+    label nm = a.n()*a.m();
 
     if (nm)
     {
         label curMinI = 0;
-        const Type* v = a.v_[0];
+        const Type* v = a[0];
 
         for (register label i=1; i<nm; i++)
         {
@@ -301,14 +301,14 @@ const Type& Foam::min(const Matrix<Form, Type>& a)
 template<class Form, class Type>
 Form Foam::operator-(const Matrix<Form, Type>& a)
 {
-    Form na(a.n_, a.m_);
+    Form na(a.n(), a.m());
 
-    if (a.n_ && a.m_)
+    if (a.n() && a.m())
     {
-        Type* nav = na.v_[0];
-        const Type* av = a.v_[0];
+        Type* nav = na[0];
+        const Type* av = a[0];
 
-        label nm = a.n_*a.m_;
+        label nm = a.n()*a.m();
         for (register label i=0; i<nm; i++)
         {
             nav[i] = -av[i];
@@ -322,33 +322,33 @@ Form Foam::operator-(const Matrix<Form, Type>& a)
 template<class Form, class Type>
 Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 {
-    if (a.n_ != b.n_)
+    if (a.n() != b.n())
     {
         FatalErrorIn
         (
             "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
         )   << "attempted add matrices with different number of rows: "
-            << a.n_ << ", " << b.n_
+            << a.n() << ", " << b.n()
             << abort(FatalError);
     }
 
-    if (a.m_ != b.m_)
+    if (a.m() != b.m())
     {
         FatalErrorIn
         (
             "Matrix<Form, Type>::operator+(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
         )   << "attempted add matrices with different number of columns: "
-            << a.m_ << ", " << b.m_
+            << a.m() << ", " << b.m()
             << abort(FatalError);
     }
 
-    Form ab(a.n_, a.m_);
+    Form ab(a.n(), a.m());
 
-    Type* abv = ab.v_[0];
-    const Type* av = a.v_[0];
-    const Type* bv = b.v_[0];
+    Type* abv = ab[0];
+    const Type* av = a[0];
+    const Type* bv = b[0];
 
-    label nm = a.n_*a.m_;
+    label nm = a.n()*a.m();
     for (register label i=0; i<nm; i++)
     {
         abv[i] = av[i] + bv[i];
@@ -361,33 +361,33 @@ Form Foam::operator+(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 template<class Form, class Type>
 Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 {
-    if (a.n_ != b.n_)
+    if (a.n() != b.n())
     {
         FatalErrorIn
         (
             "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
         )   << "attempted add matrices with different number of rows: "
-            << a.n_ << ", " << b.n_
+            << a.n() << ", " << b.n()
             << abort(FatalError);
     }
 
-    if (a.m_ != b.m_)
+    if (a.m() != b.m())
     {
         FatalErrorIn
         (
             "Matrix<Form, Type>::operator-(const Matrix<Form, Type>&, const Matrix<Form, Type>&)"
         )   << "attempted add matrices with different number of columns: "
-            << a.m_ << ", " << b.m_
+            << a.m() << ", " << b.m()
             << abort(FatalError);
     }
 
-    Form ab(a.n_, a.m_);
+    Form ab(a.n(), a.m());
 
-    Type* abv = ab.v_[0];
-    const Type* av = a.v_[0];
-    const Type* bv = b.v_[0];
+    Type* abv = ab[0];
+    const Type* av = a[0];
+    const Type* bv = b[0];
 
-    label nm = a.n_*a.m_;
+    label nm = a.n()*a.m();
     for (register label i=0; i<nm; i++)
     {
         abv[i] = av[i] - bv[i];
@@ -400,14 +400,14 @@ Form Foam::operator-(const Matrix<Form, Type>& a, const Matrix<Form, Type>& b)
 template<class Form, class Type>
 Form Foam::operator*(const scalar s, const Matrix<Form, Type>& a)
 {
-    Form sa(a.n_, a.m_);
+    Form sa(a.n(), a.m());
 
-    if (a.n_ && a.m_)
+    if (a.n() && a.m())
     {
-        Type* sav = sa.v_[0];
-        const Type* av = a.v_[0];
+        Type* sav = sa[0];
+        const Type* av = a[0];
 
-        label nm = a.n_*a.m_;
+        label nm = a.n()*a.m();
         for (register label i=0; i<nm; i++)
         {
             sav[i] = s*av[i];
diff --git a/src/OpenFOAM/matrices/Matrix/Matrix.H b/src/OpenFOAM/matrices/Matrix/Matrix.H
index fd623fdb20..b674956e2f 100644
--- a/src/OpenFOAM/matrices/Matrix/Matrix.H
+++ b/src/OpenFOAM/matrices/Matrix/Matrix.H
@@ -161,10 +161,10 @@ public:
 
     // Member operators
 
-        //- Return subscript-checked element of Matrix.
+        //- Return subscript-checked row of Matrix.
         inline Type* operator[](const label);
 
-        //- Return subscript-checked element of constant Matrix.
+        //- Return subscript-checked row of constant Matrix.
         inline const Type* operator[](const label) const;
 
         //- Assignment operator. Takes linear time.
diff --git a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
index 11be3c3b86..36e6fe6afd 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C
@@ -28,6 +28,9 @@ License
 #include "polyBoundaryMesh.H"
 #include "facePointPatch.H"
 #include "globalPointPatch.H"
+#include "PstreamBuffers.H"
+#include "lduSchedule.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -58,14 +61,46 @@ Foam::pointBoundaryMesh::pointBoundaryMesh
 
 void Foam::pointBoundaryMesh::calcGeometry()
 {
-    forAll(*this, patchi)
-    {
-        operator[](patchi).initGeometry();
-    }
+    PstreamBuffers pBufs(Pstream::defaultCommsType);
 
-    forAll(*this, patchi)
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
     {
-        operator[](patchi).calcGeometry();
+        forAll(*this, patchi)
+        {
+            operator[](patchi).initGeometry(pBufs);
+        }
+
+        pBufs.finishedSends();
+
+        forAll(*this, patchi)
+        {
+            operator[](patchi).calcGeometry(pBufs);
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule = mesh().globalData().patchSchedule();
+
+        // Dummy.
+        pBufs.finishedSends();
+
+        forAll(patchSchedule, patchEvali)
+        {
+            label patchi = patchSchedule[patchEvali].patch;
+
+            if (patchSchedule[patchEvali].init)
+            {
+                operator[](patchi).initGeometry(pBufs);
+            }
+            else
+            {
+                operator[](patchi).calcGeometry(pBufs);
+            }
+        }
     }
 }
 
@@ -97,32 +132,92 @@ Foam::pointBoundaryMesh::globalPatch() const
 
 void Foam::pointBoundaryMesh::movePoints(const pointField& p)
 {
-    pointPatchList& patches = *this;
+    PstreamBuffers pBufs(Pstream::defaultCommsType);
 
-    forAll(patches, patchi)
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
     {
-        patches[patchi].initMovePoints(p);
+        forAll(*this, patchi)
+        {
+            operator[](patchi).initMovePoints(pBufs, p);
+        }
+
+        pBufs.finishedSends();
+
+        forAll(*this, patchi)
+        {
+            operator[](patchi).movePoints(pBufs, p);
+        }
     }
-
-    forAll(patches, patchi)
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
     {
-        patches[patchi].movePoints(p);
+        const lduSchedule& patchSchedule = mesh().globalData().patchSchedule();
+
+        // Dummy.
+        pBufs.finishedSends();
+
+        forAll(patchSchedule, patchEvali)
+        {
+            label patchi = patchSchedule[patchEvali].patch;
+
+            if (patchSchedule[patchEvali].init)
+            {
+                operator[](patchi).initMovePoints(pBufs, p);
+            }
+            else
+            {
+                operator[](patchi).movePoints(pBufs, p);
+            }
+        }
     }
 }
 
 
 void Foam::pointBoundaryMesh::updateMesh()
 {
-    pointPatchList& patches = *this;
+    PstreamBuffers pBufs(Pstream::defaultCommsType);
 
-    forAll(patches, patchi)
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
     {
-        patches[patchi].initUpdateMesh();
+        forAll(*this, patchi)
+        {
+            operator[](patchi).initUpdateMesh(pBufs);
+        }
+
+        pBufs.finishedSends();
+
+        forAll(*this, patchi)
+        {
+            operator[](patchi).updateMesh(pBufs);
+        }
     }
-
-    forAll(patches, patchi)
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
     {
-        patches[patchi].updateMesh();
+        const lduSchedule& patchSchedule = mesh().globalData().patchSchedule();
+
+        // Dummy.
+        pBufs.finishedSends();
+
+        forAll(patchSchedule, patchEvali)
+        {
+            label patchi = patchSchedule[patchEvali].patch;
+
+            if (patchSchedule[patchEvali].init)
+            {
+                operator[](patchi).initUpdateMesh(pBufs);
+            }
+            else
+            {
+                operator[](patchi).updateMesh(pBufs);
+            }
+        }
     }
 }
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/basic/coupled/coupledPointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/basic/coupled/coupledPointPatch.H
index 77cbb4162f..79a2f2363c 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/basic/coupled/coupledPointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/basic/coupled/coupledPointPatch.H
@@ -66,22 +66,22 @@ protected:
     // Protected Member Functions
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry() = 0;
+        virtual void initGeometry(PstreamBuffers&) = 0;
 
         //- Calculate the patch geometry
-        virtual void calcGeometry() = 0;
+        virtual void calcGeometry(PstreamBuffers&) = 0;
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&) = 0;
+        virtual void initMovePoints(PstreamBuffers&, const pointField&) = 0;
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&) = 0;
+        virtual void movePoints(PstreamBuffers&, const pointField&) = 0;
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh() = 0;
+        virtual void initUpdateMesh(PstreamBuffers&) = 0;
 
         //- Update of the patch topology
-        virtual void updateMesh() = 0;
+        virtual void updateMesh(PstreamBuffers&) = 0;
 
 
 public:
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.C b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.C
index 81104d726e..6f65c70764 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.C
@@ -50,13 +50,13 @@ addToRunTimeSelectionTable
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
-void Foam::cyclicPointPatch::initGeometry()
+void Foam::cyclicPointPatch::initGeometry(PstreamBuffers&)
 {
     transformPairs_.setSize(0);
 }
 
 
-void Foam::cyclicPointPatch::calcGeometry()
+void Foam::cyclicPointPatch::calcGeometry(PstreamBuffers&)
 {
     const edgeList& cp = cyclicPolyPatch_.coupledPoints();
     const labelList& mp = cyclicPolyPatch_.meshPoints();
@@ -128,16 +128,20 @@ void Foam::cyclicPointPatch::calcGeometry()
             }
             else if (pointMap[cp[i][0]] == -1 && pointMap[cp[i][1]] != -1)
             {
-                FatalErrorIn("cyclicPointPatch::calcGeometry() const")
-                    << "Point " << cp[i][0] << "of point-pair " << i
+                FatalErrorIn
+                (
+                    "cyclicPointPatch::calcGeometry(PstreamBuffers&) const"
+                )   << "Point " << cp[i][0] << "of point-pair " << i
                     << " is a global point but the other point "
                     << cp[i][1] << " is not"
                     << exit(FatalError);
             }
             else if (pointMap[cp[i][0]] != -1 && pointMap[cp[i][1]] == -1)
             {
-                FatalErrorIn("cyclicPointPatch::calcGeometry() const")
-                    << "Point " << cp[i][1] << "of point-pair " << i
+                FatalErrorIn
+                (
+                    "cyclicPointPatch::calcGeometry(PstreamBuffers&) const"
+                )   << "Point " << cp[i][1] << "of point-pair " << i
                     << " is a global point but the other point "
                     << cp[i][0] << " is not"
                     << exit(FatalError);
@@ -149,25 +153,25 @@ void Foam::cyclicPointPatch::calcGeometry()
 }
 
 
-void cyclicPointPatch::initMovePoints(const pointField&)
+void cyclicPointPatch::initMovePoints(PstreamBuffers&, const pointField&)
 {}
 
 
-void cyclicPointPatch::movePoints(const pointField&)
+void cyclicPointPatch::movePoints(PstreamBuffers&, const pointField&)
 {}
 
 
-void cyclicPointPatch::initUpdateMesh()
+void cyclicPointPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    facePointPatch::initUpdateMesh();
-    cyclicPointPatch::initGeometry();
+    facePointPatch::initUpdateMesh(pBufs);
+    cyclicPointPatch::initGeometry(pBufs);
 }
 
 
-void cyclicPointPatch::updateMesh()
+void cyclicPointPatch::updateMesh(PstreamBuffers& pBufs)
 {
-    facePointPatch::updateMesh();
-    cyclicPointPatch::calcGeometry();
+    facePointPatch::updateMesh(pBufs);
+    cyclicPointPatch::calcGeometry(pBufs);
 }
 
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H
index 999c45613d..1d013f3421 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H
@@ -74,22 +74,22 @@ class cyclicPointPatch
         edgeList transformPairs_;
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry();
+        virtual void initGeometry(PstreamBuffers&);
 
         //- Calculate the patch geometry
-        virtual void calcGeometry();
+        virtual void calcGeometry(PstreamBuffers&);
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&);
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&);
+        virtual void movePoints(PstreamBuffers&, const pointField&);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh();
+        virtual void initUpdateMesh(PstreamBuffers&);
 
         //- Update of the patch topology
-        virtual void updateMesh();
+        virtual void updateMesh(PstreamBuffers&);
 
 
 public:
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
index eb99350c13..d3aae3ec87 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C
@@ -52,7 +52,7 @@ addToRunTimeSelectionTable
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
-void Foam::processorPointPatch::initGeometry()
+void Foam::processorPointPatch::initGeometry(PstreamBuffers& pBufs)
 {
     // Algorithm:
     // Depending on whether the patch is a master or a slave, get the primitive
@@ -84,16 +84,16 @@ void Foam::processorPointPatch::initGeometry()
 
     if (Pstream::parRun())
     {
-        initPatchPatchPoints();
+        initPatchPatchPoints(pBufs);
     }
 }
 
 
-void Foam::processorPointPatch::calcGeometry()
+void Foam::processorPointPatch::calcGeometry(PstreamBuffers& pBufs)
 {
     if (Pstream::parRun())
     {
-        calcPatchPatchPoints();
+        calcPatchPatchPoints(pBufs);
     }
 
     // If it is not runing parallel or there are no global points
@@ -149,11 +149,11 @@ void Foam::processorPointPatch::calcGeometry()
 }
 
 
-void processorPointPatch::initPatchPatchPoints()
+void processorPointPatch::initPatchPatchPoints(PstreamBuffers& pBufs)
 {
     if (debug)
     {
-        Info<< "processorPointPatch::calcPatchPatchPoints() : "
+        Info<< "processorPointPatch::initPatchPatchPoints(PstreamBuffers&) : "
             << "constructing patch-patch points"
             << endl;
     }
@@ -229,7 +229,7 @@ void processorPointPatch::initPatchPatchPoints()
 
     // Send the patchPatchPoints to the neighbouring processor
 
-    OPstream toNeighbProc(Pstream::blocking, neighbProcNo());
+    UOPstream toNeighbProc(neighbProcNo(), pBufs);
 
     toNeighbProc
         << ppmp.size()              // number of points for checking
@@ -238,17 +238,17 @@ void processorPointPatch::initPatchPatchPoints()
 
     if (debug)
     {
-        Info<< "processorPointPatch::calcPatchPatchPoints() : "
+        Info<< "processorPointPatch::initPatchPatchPoints() : "
             << "constructed patch-patch points"
             << endl;
     }
 }
 
 
-void Foam::processorPointPatch::calcPatchPatchPoints()
+void Foam::processorPointPatch::calcPatchPatchPoints(PstreamBuffers& pBufs)
 {
     // Get the patchPatchPoints from the neighbouring processor
-    IPstream fromNeighbProc(Pstream::blocking, neighbProcNo());
+    UIPstream fromNeighbProc(neighbProcNo(), pBufs);
 
     label nbrNPoints(readLabel(fromNeighbProc));
     labelListList patchPatchPoints(fromNeighbProc);
@@ -265,7 +265,7 @@ void Foam::processorPointPatch::calcPatchPatchPoints()
     // separate.
     if (nbrNPoints != ppmp.size())
     {
-        WarningIn("processorPointPatch::calcPatchPatchPoints()")
+        WarningIn("processorPointPatch::calcPatchPatchPoints(PstreamBuffers&)")
             << "Processor patch " << name()
             << " has " << ppmp.size() << " points; coupled patch has "
             << nbrNPoints << " points." << endl
@@ -352,25 +352,25 @@ void Foam::processorPointPatch::calcPatchPatchPoints()
 }
 
 
-void processorPointPatch::initMovePoints(const pointField&)
+void processorPointPatch::initMovePoints(PstreamBuffers&, const pointField&)
 {}
 
 
-void processorPointPatch::movePoints(const pointField&)
+void processorPointPatch::movePoints(PstreamBuffers&, const pointField&)
 {}
 
 
-void processorPointPatch::initUpdateMesh()
+void processorPointPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    facePointPatch::initUpdateMesh();
-    processorPointPatch::initGeometry();
+    facePointPatch::initUpdateMesh(pBufs);
+    processorPointPatch::initGeometry(pBufs);
 }
 
 
-void processorPointPatch::updateMesh()
+void processorPointPatch::updateMesh(PstreamBuffers& pBufs)
 {
-    facePointPatch::updateMesh();
-    processorPointPatch::calcGeometry();
+    facePointPatch::updateMesh(pBufs);
+    processorPointPatch::calcGeometry(pBufs);
 }
 
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.H
index f17a00ca83..373416ac0a 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.H
@@ -68,30 +68,30 @@ class processorPointPatch
     // Private Member Functions
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry();
+        virtual void initGeometry(PstreamBuffers&);
 
         //- Calculate the patch geometry
-        virtual void calcGeometry();
+        virtual void calcGeometry(PstreamBuffers&);
 
         //- Initialise the points on this patch which are should also be 
         //  on a neighbouring patch but are not part of faces of that patch
-        void initPatchPatchPoints();
+        void initPatchPatchPoints(PstreamBuffers&);
 
         //- Calculate the points on this patch which are should also be 
         //  on a neighbouring patch but are not part of faces of that patch
-        void calcPatchPatchPoints();
+        void calcPatchPatchPoints(PstreamBuffers&);
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&);
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&);
+        virtual void movePoints(PstreamBuffers&, const pointField&);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh();
+        virtual void initUpdateMesh(PstreamBuffers&);
 
         //- Update of the patch topology
-        virtual void updateMesh();
+        virtual void updateMesh(PstreamBuffers&);
 
 
         //- Disallow default construct as copy
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H
index dc2437877d..a77af7af30 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H
@@ -86,7 +86,7 @@ protected:
         // Construction of demand-driven data
 
             //- Calculate mesh points
-            virtual void calcGeometry() = 0;
+            virtual void calcGeometry(PstreamBuffers&) = 0;
 
 
 public:
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/global/globalPointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/global/globalPointPatch.H
index 990ab0de9a..ba89ca950c 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/global/globalPointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/derived/global/globalPointPatch.H
@@ -67,27 +67,27 @@ class globalPointPatch
     // Protected Member Functions
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry()
+        virtual void initGeometry(PstreamBuffers&)
         {}
 
         //- Calculate the patch geometry
-        virtual void calcGeometry()
+        virtual void calcGeometry(PstreamBuffers&)
         {}
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&)
+        virtual void initMovePoints(PstreamBuffers&, const pointField&)
         {}
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&)
+        virtual void movePoints(PstreamBuffers&, const pointField&)
         {}
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh()
+        virtual void initUpdateMesh(PstreamBuffers&)
         {}
 
         //- Update of the patch topology
-        virtual void updateMesh()
+        virtual void updateMesh(PstreamBuffers&)
         {}
 
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.C b/src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.C
index 25c91dd398..fb15ac4b96 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.C
@@ -51,7 +51,7 @@ addToRunTimeSelectionTable
 
 // * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
 
-void facePointPatch::initGeometry()
+void facePointPatch::initGeometry(PstreamBuffers&)
 {
     meshPoints_.setSize(0);
     localPoints_.setSize(0);
@@ -59,25 +59,25 @@ void facePointPatch::initGeometry()
 }
 
 
-void facePointPatch::calcGeometry()
+void facePointPatch::calcGeometry(PstreamBuffers&)
 {}
 
 
-void facePointPatch::initMovePoints(const pointField&)
+void facePointPatch::initMovePoints(PstreamBuffers&, const pointField&)
 {}
 
 
-void facePointPatch::movePoints(const pointField&)
+void facePointPatch::movePoints(PstreamBuffers&, const pointField&)
 {}
 
 
-void facePointPatch::initUpdateMesh()
+void facePointPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    facePointPatch::initGeometry();
+    facePointPatch::initGeometry(pBufs);
 }
 
 
-void facePointPatch::updateMesh()
+void facePointPatch::updateMesh(PstreamBuffers&)
 {}
 
 
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.H
index 7a75b4017f..0094d4ce1a 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.H
@@ -76,22 +76,22 @@ protected:
     // Protected Member Functions
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry();
+        virtual void initGeometry(PstreamBuffers&);
 
         //- Calculate the patch geometry
-        virtual void calcGeometry();
+        virtual void calcGeometry(PstreamBuffers&);
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&);
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&);
+        virtual void movePoints(PstreamBuffers&, const pointField&);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh();
+        virtual void initUpdateMesh(PstreamBuffers&);
 
         //- Update of the patch topology
-        virtual void updateMesh();
+        virtual void updateMesh(PstreamBuffers&);
 
 
 private:
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/pointPatch/pointPatch.H b/src/OpenFOAM/meshes/pointMesh/pointPatches/pointPatch/pointPatch.H
index bdf67c879b..637ac659dd 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/pointPatch/pointPatch.H
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/pointPatch/pointPatch.H
@@ -49,6 +49,7 @@ namespace Foam
 
 class pointBoundaryMesh;
 class pointConstraint;
+class PstreamBuffers;
 
 /*---------------------------------------------------------------------------*\
                       Class pointPatch Declaration
@@ -79,27 +80,27 @@ protected:
         friend class pointBoundaryMesh;
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry()
+        virtual void initGeometry(PstreamBuffers&)
         {}
 
         //- Calculate the patch geometry
-        virtual void calcGeometry()
+        virtual void calcGeometry(PstreamBuffers&)
         {}
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&)
+        virtual void initMovePoints(PstreamBuffers&, const pointField&)
         {}
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&)
+        virtual void movePoints(PstreamBuffers&, const pointField&)
         {}
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh()
+        virtual void initUpdateMesh(PstreamBuffers&)
         {}
 
         //- Update of the patch topology
-        virtual void updateMesh()
+        virtual void updateMesh(PstreamBuffers&)
         {}
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index eff7af933c..a111e0435d 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -29,6 +29,9 @@ License
 #include "primitiveMesh.H"
 #include "processorPolyPatch.H"
 #include "stringListOps.H"
+#include "PstreamBuffers.H"
+#include "lduSchedule.H"
+#include "globalMeshData.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -144,14 +147,46 @@ void Foam::polyBoundaryMesh::clearAddressing()
 
 void Foam::polyBoundaryMesh::calcGeometry()
 {
-    forAll(*this, patchi)
-    {
-        operator[](patchi).initGeometry();
-    }
+    PstreamBuffers pBufs(Pstream::defaultCommsType);
 
-    forAll(*this, patchi)
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
     {
-        operator[](patchi).calcGeometry();
+        forAll(*this, patchi)
+        {
+            operator[](patchi).initGeometry(pBufs);
+        }
+
+        pBufs.finishedSends();
+
+        forAll(*this, patchi)
+        {
+            operator[](patchi).calcGeometry(pBufs);
+        }
+    }
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
+    {
+        const lduSchedule& patchSchedule = mesh().globalData().patchSchedule();
+
+        // Dummy.
+        pBufs.finishedSends();
+
+        forAll(patchSchedule, patchEvali)
+        {
+            label patchi = patchSchedule[patchEvali].patch;
+
+            if (patchSchedule[patchEvali].init)
+            {
+                operator[](patchi).initGeometry(pBufs);
+            }
+            else
+            {
+                operator[](patchi).calcGeometry(pBufs);
+            }
+        }
     }
 }
 
@@ -573,16 +608,46 @@ bool Foam::polyBoundaryMesh::checkDefinition(const bool report) const
 
 void Foam::polyBoundaryMesh::movePoints(const pointField& p)
 {
-    polyPatchList& patches = *this;
+    PstreamBuffers pBufs(Pstream::defaultCommsType);
 
-    forAll(patches, patchi)
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
     {
-        patches[patchi].initMovePoints(p);
+        forAll(*this, patchi)
+        {
+            operator[](patchi).initMovePoints(pBufs, p);
+        }
+
+        pBufs.finishedSends();
+
+        forAll(*this, patchi)
+        {
+            operator[](patchi).movePoints(pBufs, p);
+        }
     }
-
-    forAll(patches, patchi)
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
     {
-        patches[patchi].movePoints(p);
+        const lduSchedule& patchSchedule = mesh().globalData().patchSchedule();
+
+        // Dummy.
+        pBufs.finishedSends();
+
+        forAll(patchSchedule, patchEvali)
+        {
+            label patchi = patchSchedule[patchEvali].patch;
+
+            if (patchSchedule[patchEvali].init)
+            {
+                operator[](patchi).initMovePoints(pBufs, p);
+            }
+            else
+            {
+                operator[](patchi).movePoints(pBufs, p);
+            }
+        }
     }
 }
 
@@ -591,16 +656,46 @@ void Foam::polyBoundaryMesh::updateMesh()
 {
     deleteDemandDrivenData(neighbourEdgesPtr_);
 
-    polyPatchList& patches = *this;
+    PstreamBuffers pBufs(Pstream::defaultCommsType);
 
-    forAll(patches, patchi)
+    if
+    (
+        Pstream::defaultCommsType == Pstream::blocking
+     || Pstream::defaultCommsType == Pstream::nonBlocking
+    )
     {
-        patches[patchi].initUpdateMesh();
+        forAll(*this, patchi)
+        {
+            operator[](patchi).initUpdateMesh(pBufs);
+        }
+
+        pBufs.finishedSends();
+
+        forAll(*this, patchi)
+        {
+            operator[](patchi).updateMesh(pBufs);
+        }
     }
-
-    forAll(patches, patchi)
+    else if (Pstream::defaultCommsType == Pstream::scheduled)
     {
-        patches[patchi].updateMesh();
+        const lduSchedule& patchSchedule = mesh().globalData().patchSchedule();
+
+        // Dummy.
+        pBufs.finishedSends();
+
+        forAll(patchSchedule, patchEvali)
+        {
+            label patchi = patchSchedule[patchEvali].patch;
+
+            if (patchSchedule[patchEvali].init)
+            {
+                operator[](patchi).initUpdateMesh(pBufs);
+            }
+            else
+            {
+                operator[](patchi).updateMesh(pBufs);
+            }
+        }
     }
 }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index e3277defcc..fdea473ed3 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
@@ -50,8 +50,6 @@ class polyMesh;
 
 // Forward declaration of friend functions and operators
 
-class polyBoundaryMesh;
-
 Ostream& operator<<(Ostream&, const polyBoundaryMesh&);
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
index f823148d3a..e1afad60a8 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
@@ -90,22 +90,22 @@ protected:
         ) const;
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry() = 0;
+        virtual void initGeometry(PstreamBuffers&) = 0;
 
         //- Calculate the patch geometry
-        virtual void calcGeometry() = 0;
+        virtual void calcGeometry(PstreamBuffers&) = 0;
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&) = 0;
+        virtual void initMovePoints(PstreamBuffers&, const pointField&) = 0;
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&) = 0;
+        virtual void movePoints(PstreamBuffers&, const pointField&) = 0;
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh() = 0;
+        virtual void initUpdateMesh(PstreamBuffers&) = 0;
 
         //- Update of the patch topology
-        virtual void updateMesh() = 0;
+        virtual void updateMesh(PstreamBuffers&) = 0;
 
 
         //- Write point in OBJ format
@@ -283,7 +283,11 @@ public:
 
         //- Initialize ordering for primitivePatch. Does not
         //  refer to *this (except for name() and type() etc.)
-        virtual void initOrder(const primitivePatch&) const = 0;
+        virtual void initOrder
+        (
+            PstreamBuffers&,
+            const primitivePatch&
+        ) const = 0;
 
         //- Return new ordering for primitivePatch.
         //  Ordering is -faceMap: for every face
@@ -292,6 +296,7 @@ public:
         //  (faceMap is identity, rotation is 0), true otherwise.
         virtual bool order
         (
+            PstreamBuffers&,
             const primitivePatch&,
             labelList& faceMap,
             labelList& rotation
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
index 8c6b202079..fef6c98ec0 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.C
@@ -516,7 +516,8 @@ void Foam::cyclicPolyPatch::getCentresAndAnchors
         //    if (debug)
         //    {
         //        Pout<< "cyclicPolyPatch::getCentresAndAnchors :"
-        //            << "Specified translation : " << separationVector_ << endl;
+        //            << "Specified translation : " << separationVector_
+        //            << endl;
         //    }
         //
         //    half0Ctrs += separationVector_;
@@ -858,36 +859,44 @@ Foam::cyclicPolyPatch::~cyclicPolyPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::cyclicPolyPatch::initGeometry()
+void Foam::cyclicPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
-    polyPatch::initGeometry();
+    polyPatch::initGeometry(pBufs);
 }
 
-void Foam::cyclicPolyPatch::calcGeometry()
+void Foam::cyclicPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 {
-    polyPatch::calcGeometry();
+    polyPatch::calcGeometry(pBufs);
     calcTransforms();
 }
 
-void Foam::cyclicPolyPatch::initMovePoints(const pointField& p)
+void Foam::cyclicPolyPatch::initMovePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
 {
-    polyPatch::initMovePoints(p);
+    polyPatch::initMovePoints(pBufs, p);
 }
 
-void Foam::cyclicPolyPatch::movePoints(const pointField& p)
+void Foam::cyclicPolyPatch::movePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
 {
-    polyPatch::movePoints(p);
+    polyPatch::movePoints(pBufs, p);
     calcTransforms();
 }
 
-void Foam::cyclicPolyPatch::initUpdateMesh()
+void Foam::cyclicPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    polyPatch::initUpdateMesh();
+    polyPatch::initUpdateMesh(pBufs);
 }
 
-void Foam::cyclicPolyPatch::updateMesh()
+void Foam::cyclicPolyPatch::updateMesh(PstreamBuffers& pBufs)
 {
-    polyPatch::updateMesh();
+    polyPatch::updateMesh(pBufs);
     deleteDemandDrivenData(coupledPointsPtr_);
     deleteDemandDrivenData(coupledEdgesPtr_);
 }
@@ -1105,7 +1114,11 @@ const Foam::edgeList& Foam::cyclicPolyPatch::coupledEdges() const
 }
 
 
-void Foam::cyclicPolyPatch::initOrder(const primitivePatch& pp) const
+void Foam::cyclicPolyPatch::initOrder
+(
+    PstreamBuffers&,
+    const primitivePatch& pp
+) const
 {}
 
 
@@ -1115,6 +1128,7 @@ void Foam::cyclicPolyPatch::initOrder(const primitivePatch& pp) const
 //  is identity, rotation is 0)
 bool Foam::cyclicPolyPatch::order
 (
+    PstreamBuffers& pBufs,
     const primitivePatch& pp,
     labelList& faceMap,
     labelList& rotation
@@ -1300,7 +1314,8 @@ bool Foam::cyclicPolyPatch::order
                 << endl;
 
             // Recalculate untransformed face centres
-            //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+            //pointField rawHalf0Ctrs =
+            //    calcFaceCentres(half0Faces, pp.points());
             label vertI = 0;
 
             forAll(half1Ctrs, i)
@@ -1413,7 +1428,8 @@ bool Foam::cyclicPolyPatch::order
                     << endl;
 
                 // Recalculate untransformed face centres
-                //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+                //pointField rawHalf0Ctrs =
+                //    calcFaceCentres(half0Faces, pp.points());
                 label vertI = 0;
 
                 forAll(half1Ctrs, i)
@@ -1499,7 +1515,8 @@ bool Foam::cyclicPolyPatch::order
                 << endl;
 
             // Recalculate untransformed face centres
-            //pointField rawHalf0Ctrs = calcFaceCentres(half0Faces, pp.points());
+            //pointField rawHalf0Ctrs =
+            //    calcFaceCentres(half0Faces, pp.points());
             label vertI = 0;
 
             forAll(half1Ctrs, i)
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
index ea67348891..4f81975e60 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/cyclic/cyclicPolyPatch.H
@@ -176,22 +176,22 @@ protected:
     // Protected Member functions
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry();
+        virtual void initGeometry(PstreamBuffers&);
 
         //- Calculate the patch geometry
-        virtual void calcGeometry();
+        virtual void calcGeometry(PstreamBuffers&);
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&);
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&);
+        virtual void movePoints(PstreamBuffers&, const pointField&);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh();
+        virtual void initUpdateMesh(PstreamBuffers&);
 
         //- Update of the patch topology
-        virtual void updateMesh();
+        virtual void updateMesh(PstreamBuffers&);
 
 public:
 
@@ -391,7 +391,7 @@ public:
 
         //- Initialize ordering for primitivePatch. Does not
         //  refer to *this (except for name() and type() etc.)
-        virtual void initOrder(const primitivePatch&) const;
+        virtual void initOrder(PstreamBuffers&, const primitivePatch&) const;
 
         //- Return new ordering for primitivePatch.
         //  Ordering is -faceMap: for every face
@@ -400,6 +400,7 @@ public:
         //  (faceMap is identity, rotation is 0), true otherwise.
         virtual bool order
         (
+            PstreamBuffers&,
             const primitivePatch&,
             labelList& faceMap,
             labelList& rotation
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
index 6e4df9f4e5..05969ee593 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.C
@@ -34,6 +34,7 @@ License
 #include "polyMesh.H"
 #include "Time.H"
 #include "transformList.H"
+#include "PstreamBuffers.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -155,16 +156,11 @@ Foam::processorPolyPatch::~processorPolyPatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::processorPolyPatch::initGeometry()
+void Foam::processorPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
     if (Pstream::parRun())
     {
-        OPstream toNeighbProc
-        (
-            Pstream::blocking,
-            neighbProcNo(),
-            3*(sizeof(label) + size()*sizeof(vector) + sizeof(scalar))
-        );
+        UOPstream toNeighbProc(neighbProcNo(), pBufs);
 
         toNeighbProc
             << faceCentres()
@@ -174,17 +170,13 @@ void Foam::processorPolyPatch::initGeometry()
 }
 
 
-void Foam::processorPolyPatch::calcGeometry()
+void Foam::processorPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 {
     if (Pstream::parRun())
     {
         {
-            IPstream fromNeighbProc
-            (
-                Pstream::blocking,
-                neighbProcNo(),
-                3*(sizeof(label) + size()*sizeof(vector) + sizeof(scalar))
-            );
+            UIPstream fromNeighbProc(neighbProcNo(), pBufs);
+
             fromNeighbProc
                 >> neighbFaceCentres_
                 >> neighbFaceAreas_
@@ -251,22 +243,30 @@ void Foam::processorPolyPatch::calcGeometry()
 }
 
 
-void Foam::processorPolyPatch::initMovePoints(const pointField& p)
+void Foam::processorPolyPatch::initMovePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
 {
-    polyPatch::movePoints(p);
-    processorPolyPatch::initGeometry();
+    polyPatch::movePoints(pBufs, p);
+    processorPolyPatch::initGeometry(pBufs);
 }
 
 
-void Foam::processorPolyPatch::movePoints(const pointField&)
+void Foam::processorPolyPatch::movePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField&
+)
 {
-    processorPolyPatch::calcGeometry();
+    processorPolyPatch::calcGeometry(pBufs);
 }
 
 
-void Foam::processorPolyPatch::initUpdateMesh()
+void Foam::processorPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    polyPatch::initUpdateMesh();
+    polyPatch::initUpdateMesh(pBufs);
 
     deleteDemandDrivenData(neighbPointsPtr_);
     deleteDemandDrivenData(neighbEdgesPtr_);
@@ -303,14 +303,7 @@ void Foam::processorPolyPatch::initUpdateMesh()
             edgeIndex[patchEdgeI] = findIndex(fEdges, patchEdgeI);
         }
 
-        OPstream toNeighbProc
-        (
-            Pstream::blocking,
-            neighbProcNo(),
-            8*sizeof(label)             // four headers of labelList
-          + 2*nPoints()*sizeof(label)   // two point-based labellists
-          + 2*nEdges()*sizeof(label)    // two edge-based labelLists
-        );
+        UOPstream toNeighbProc(neighbProcNo(), pBufs);
 
         toNeighbProc
             << pointFace
@@ -321,10 +314,10 @@ void Foam::processorPolyPatch::initUpdateMesh()
 }
 
 
-void Foam::processorPolyPatch::updateMesh()
+void Foam::processorPolyPatch::updateMesh(PstreamBuffers& pBufs)
 {
     // For completeness
-    polyPatch::updateMesh();
+    polyPatch::updateMesh(pBufs);
 
     if (Pstream::parRun())
     {
@@ -336,7 +329,7 @@ void Foam::processorPolyPatch::updateMesh()
         {
             // Note cannot predict exact size since opposite nPoints might
             // be different from one over here.
-            IPstream fromNeighbProc(Pstream::blocking, neighbProcNo());
+            UIPstream fromNeighbProc(neighbProcNo(), pBufs);
 
             fromNeighbProc
                 >> nbrPointFace
@@ -446,7 +439,11 @@ const Foam::labelList& Foam::processorPolyPatch::neighbEdges() const
 }
 
 
-void Foam::processorPolyPatch::initOrder(const primitivePatch& pp) const
+void Foam::processorPolyPatch::initOrder
+(
+    PstreamBuffers& pBufs,
+    const primitivePatch& pp
+) const
 {
     if (!Pstream::parRun())
     {
@@ -491,7 +488,7 @@ void Foam::processorPolyPatch::initOrder(const primitivePatch& pp) const
         pointField anchors(getAnchorPoints(pp, pp.points()));
 
         // Now send all info over to the neighbour
-        OPstream toNeighbour(Pstream::blocking, neighbProcNo());
+        UOPstream toNeighbour(neighbProcNo(), pBufs);
         toNeighbour << ctrs << anchors;
     }
 }
@@ -503,6 +500,7 @@ void Foam::processorPolyPatch::initOrder(const primitivePatch& pp) const
 // is identity, rotation is 0)
 bool Foam::processorPolyPatch::order
 (
+    PstreamBuffers& pBufs,
     const primitivePatch& pp,
     labelList& faceMap,
     labelList& rotation
@@ -539,7 +537,7 @@ bool Foam::processorPolyPatch::order
 
         // Receive data from neighbour
         {
-            IPstream fromNeighbour(Pstream::blocking, neighbProcNo());
+            UIPstream fromNeighbour(neighbProcNo(), pBufs);
             fromNeighbour >> masterCtrs >> masterAnchors;
         }
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
index a8948fee43..1955892727 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/constraint/processor/processorPolyPatch.H
@@ -28,11 +28,9 @@ Class
 Description
     Neighbour processor patch.
 
-    Note: morph patch face ordering comes geometric or topological.
-    Geometric: no cyclics allowed (assumes faces coincident)
-    Topological: needs unmodified faces on both sides to correspond. Also
-    needs at least one per connected patch area (so all patch faces can be
-    visited from an unmodified face)
+    Note: morph patch face ordering tries to do a geometric ordering.
+    (assumes faces coincident) Hence will have problems when cyclics
+    are present. 
 
 SourceFiles
     processorPolyPatch.C
@@ -97,22 +95,22 @@ protected:
     // Protected Member functions
 
         //- Initialise the calculation of the patch geometry
-        void initGeometry();
+        void initGeometry(PstreamBuffers&);
 
         //- Calculate the patch geometry
-        void calcGeometry();
+        void calcGeometry(PstreamBuffers&);
 
         //- Initialise the patches for moving points
-        void initMovePoints(const pointField&);
+        void initMovePoints(PstreamBuffers&, const pointField&);
 
         //- Correct patches after moving points
-        void movePoints(const pointField&);
+        void movePoints(PstreamBuffers&, const pointField&);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh();
+        virtual void initUpdateMesh(PstreamBuffers&);
 
         //- Update of the patch topology
-        virtual void updateMesh();
+        virtual void updateMesh(PstreamBuffers&);
 
 
 public:
@@ -283,7 +281,7 @@ public:
 
         //- Initialize ordering for primitivePatch. Does not
         //  refer to *this (except for name() and type() etc.)
-        virtual void initOrder(const primitivePatch&) const;
+        virtual void initOrder(PstreamBuffers&, const primitivePatch&) const;
 
         //- Return new ordering for primitivePatch.
         //  Ordering is -faceMap: for every face
@@ -292,6 +290,7 @@ public:
         //  (faceMap is identity, rotation is 0), true otherwise.
         virtual bool order
         (
+            PstreamBuffers&,
             const primitivePatch&,
             labelList& faceMap,
             labelList& rotation
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
index be0aa378a6..0a6bde6b4f 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
@@ -55,12 +55,12 @@ namespace Foam
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-void Foam::polyPatch::movePoints(const pointField& p)
+void Foam::polyPatch::movePoints(PstreamBuffers&, const pointField& p)
 {
     primitivePatch::movePoints(p);
 }
 
-void Foam::polyPatch::updateMesh()
+void Foam::polyPatch::updateMesh(PstreamBuffers&)
 {
     clearAddressing();
 }
@@ -334,12 +334,13 @@ void Foam::polyPatch::write(Ostream& os) const
 }
 
 
-void Foam::polyPatch::initOrder(const primitivePatch&) const
+void Foam::polyPatch::initOrder(PstreamBuffers&, const primitivePatch&) const
 {}
 
 
 bool Foam::polyPatch::order
 (
+    PstreamBuffers&,
     const primitivePatch&,
     labelList& faceMap,
     labelList& rotation
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
index f8c8fa72e7..5694e53823 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
@@ -56,6 +56,7 @@ namespace Foam
 
 class polyBoundaryMesh;
 class polyPatch;
+class PstreamBuffers;
 
 Ostream& operator<<(Ostream&, const polyPatch&);
 
@@ -101,26 +102,26 @@ protected:
         friend class polyBoundaryMesh;
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry()
+        virtual void initGeometry(PstreamBuffers&)
         {}
 
         //- Calculate the patch geometry
-        virtual void calcGeometry()
+        virtual void calcGeometry(PstreamBuffers&)
         {}
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&)
+        virtual void initMovePoints(PstreamBuffers&, const pointField&)
         {}
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField& p);
+        virtual void movePoints(PstreamBuffers&, const pointField& p);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh()
+        virtual void initUpdateMesh(PstreamBuffers&)
         {}
 
         //- Update of the patch topology
-        virtual void updateMesh();
+        virtual void updateMesh(PstreamBuffers&);
 
 
 public:
@@ -358,7 +359,7 @@ public:
 
         //- Initialize ordering for primitivePatch. Does not
         //  refer to *this (except for name() and type() etc.)
-        virtual void initOrder(const primitivePatch&) const;
+        virtual void initOrder(PstreamBuffers&, const primitivePatch&) const;
 
         //- Return new ordering for primitivePatch.
         //  Ordering is -faceMap: for every face
@@ -367,6 +368,7 @@ public:
         //  (faceMap is identity, rotation is 0), true otherwise.
         virtual bool order
         (
+            PstreamBuffers&,
             const primitivePatch&,
             labelList& faceMap,
             labelList& rotation
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
index 6be9ef348b..a86dd36010 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.C
@@ -28,7 +28,6 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-// Construct from components
 template
 <
     class Face,
@@ -36,7 +35,6 @@ template
     class PointField,
     class PointType
 >
-
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 PrimitivePatch
 (
@@ -60,12 +58,12 @@ PrimitivePatch
     edgeLoopsPtr_(NULL),
     localPointsPtr_(NULL),
     localPointOrderPtr_(NULL),
+    faceCentresPtr_(NULL),
     faceNormalsPtr_(NULL),
     pointNormalsPtr_(NULL)
 {}
 
 
-// Construct from components
 template
 <
     class Face,
@@ -73,7 +71,6 @@ template
     class PointField,
     class PointType
 >
-
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 PrimitivePatch
 (
@@ -98,12 +95,12 @@ PrimitivePatch
     edgeLoopsPtr_(NULL),
     localPointsPtr_(NULL),
     localPointOrderPtr_(NULL),
+    faceCentresPtr_(NULL),
     faceNormalsPtr_(NULL),
     pointNormalsPtr_(NULL)
 {}
 
 
-// Construct as copy
 template
 <
     class Face,
@@ -111,7 +108,6 @@ template
     class PointField,
     class PointType
 >
-
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 PrimitivePatch
 (
@@ -135,6 +131,7 @@ PrimitivePatch
     edgeLoopsPtr_(NULL),
     localPointsPtr_(NULL),
     localPointOrderPtr_(NULL),
+    faceCentresPtr_(NULL),
     faceNormalsPtr_(NULL),
     pointNormalsPtr_(NULL)
 {}
@@ -149,8 +146,8 @@ template
     class PointField,
     class PointType
 >
-
-Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::~PrimitivePatch()
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
+~PrimitivePatch()
 {
     clearOut();
 }
@@ -158,7 +155,6 @@ Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::~PrimitivePatch()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Correct patch after moving points
 template
 <
     class Face,
@@ -166,7 +162,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 movePoints
@@ -193,7 +188,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::edgeList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 edges() const
@@ -214,7 +208,6 @@ template
     class PointField,
     class PointType
 >
-
 Foam::label
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 nInternalEdges() const
@@ -235,7 +228,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 boundaryPoints() const
@@ -256,7 +248,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelListList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 faceFaces() const
@@ -277,7 +268,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelListList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 edgeFaces() const
@@ -298,7 +288,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelListList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 faceEdges() const
@@ -319,7 +308,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelListList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 pointEdges() const
@@ -340,7 +328,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelListList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 pointFaces() const
@@ -361,7 +348,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::List<Face>&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 localFaces() const
@@ -382,7 +368,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 meshPoints() const
@@ -403,7 +388,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::Map<Foam::label>&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 meshPointMap() const
@@ -424,7 +408,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::Field<PointType>&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 localPoints() const
@@ -445,7 +428,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 localPointOrder() const
@@ -466,7 +448,6 @@ template
     class PointField,
     class PointType
 >
-
 Foam::label
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 whichPoint
@@ -495,7 +476,26 @@ template
     class PointField,
     class PointType
 >
+const Foam::Field<PointType>&
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
+faceCentres() const
+{
+    if (!faceCentresPtr_)
+    {
+        calcFaceCentres();
+    }
 
+    return *faceCentresPtr_;
+}
+
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
 const Foam::Field<PointType>&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 faceNormals() const
@@ -516,7 +516,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::Field<PointType>&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 pointNormals() const
@@ -539,7 +538,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 operator=
@@ -552,6 +550,7 @@ operator=
     FaceList<Face>::operator=(pp);
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "PrimitivePatchAddressing.C"
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
index c741b0839e..dbd8a9ffb5 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
@@ -166,6 +166,9 @@ private:
         //- Local point order for most efficient search
         mutable labelList* localPointOrderPtr_;
 
+        //- Face centres
+        mutable Field<PointType>* faceCentresPtr_;
+
         //- Face unit normals
         mutable Field<PointType>* faceNormalsPtr_;
 
@@ -205,6 +208,9 @@ private:
         //- Calculate local point order
         void calcLocalPointOrder() const;
 
+        //- Calculate face centres
+        void calcFaceCentres() const;
+
         //- Calculate unit face normals
         void calcFaceNormals() const;
 
@@ -370,6 +376,9 @@ public:
                 const labelListList& pointEdges
             ) const;
 
+            //- Return face centres for patch
+            const Field<PointType>& faceCentres() const;
+
             //- Return face normals for patch
             const Field<PointType>& faceNormals() const;
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
index 14d91344d1..5deecba26c 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
@@ -48,7 +48,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcAddressing() const
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
index fbd309d2cf..379445b1c4 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
@@ -22,8 +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
-
 \*---------------------------------------------------------------------------*/
 
 #include "PrimitivePatch.H"
@@ -39,7 +37,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcBdryPoints() const
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchCheck.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchCheck.C
index 172bc9db7e..11fe2235ba 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchCheck.C
@@ -41,7 +41,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 visitPointRegion
@@ -120,7 +119,6 @@ template
     class PointField,
     class PointType
 >
-
 typename Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::surfaceTopo
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 surfaceType() const
@@ -174,7 +172,6 @@ template
     class PointField,
     class PointType
 >
-
 bool
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 checkTopology
@@ -245,7 +242,6 @@ template
     class PointField,
     class PointType
 >
-
 bool
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 checkPointManifold
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C
index 160bf5546f..c85cf14789 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchClear.C
@@ -22,8 +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
-
 \*---------------------------------------------------------------------------*/
 
 #include "PrimitivePatch.H"
@@ -39,7 +37,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 clearGeom()
@@ -52,6 +49,7 @@ clearGeom()
     }
 
     deleteDemandDrivenData(localPointsPtr_);
+    deleteDemandDrivenData(faceCentresPtr_);
     deleteDemandDrivenData(faceNormalsPtr_);
     deleteDemandDrivenData(pointNormalsPtr_);
 }
@@ -64,7 +62,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 clearTopology()
@@ -106,7 +103,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 clearPatchMeshAddr()
@@ -132,7 +128,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 clearOut()
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchEdgeLoops.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchEdgeLoops.C
index 3b7d87df90..80e36f166d 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchEdgeLoops.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchEdgeLoops.C
@@ -31,7 +31,6 @@ Description
 #include "PrimitivePatch.H"
 
 
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template
@@ -41,7 +40,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcEdgeLoops() const
@@ -174,7 +172,6 @@ template
     class PointField,
     class PointType
 >
-
 const Foam::labelListList&
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 edgeLoops() const
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchLocalPointOrder.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchLocalPointOrder.C
index b4359d83eb..7a20d8a81b 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchLocalPointOrder.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchLocalPointOrder.C
@@ -39,14 +39,13 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcLocalPointOrder() const
 {
     // Note: Cannot use bandCompressing as point-point addressing does
     // not exist and is not considered generally useful.
-    //  
+    //
 
     if (debug)
     {
@@ -137,4 +136,5 @@ calcLocalPointOrder() const
     }
 }
 
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
index 30cffb6310..cd05f3f0e1 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshData.C
@@ -27,8 +27,6 @@ License
 #include "PrimitivePatch.H"
 #include "Map.H"
 
-
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template
@@ -38,7 +36,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcMeshData() const
@@ -76,13 +73,13 @@ calcMeshData() const
 
     ////- 1.5 code:
     //// if the point is used, set the mark to 1
-    //forAll (*this, faceI)
+    //forAll(*this, facei)
     //{
-    //    const Face& curPoints = this->operator[](faceI);
+    //    const Face& curPoints = this->operator[](facei);
     //
-    //    forAll (curPoints, pointI)
+    //    forAll(curPoints, pointi)
     //    {
-    //        markedPoints.insert(curPoints[pointI], -1);
+    //        markedPoints.insert(curPoints[pointi], -1);
     //    }
     //}
     //
@@ -95,22 +92,22 @@ calcMeshData() const
     //sort(pointPatch);
     //
     //// For every point in map give it its label in mesh points
-    //forAll (pointPatch, pointI)
+    //forAll(pointPatch, pointi)
     //{
-    //    markedPoints.find(pointPatch[pointI])() = pointI;
+    //    markedPoints.find(pointPatch[pointi])() = pointi;
     //}
 
     //- Unsorted version:
     DynamicList<label> meshPoints(2*this->size());
-    forAll (*this, faceI)
+    forAll(*this, facei)
     {
-        const Face& curPoints = this->operator[](faceI);
+        const Face& curPoints = this->operator[](facei);
 
-        forAll (curPoints, pointI)
+        forAll(curPoints, pointi)
         {
-            if (markedPoints.insert(curPoints[pointI], meshPoints.size()))
+            if (markedPoints.insert(curPoints[pointi], meshPoints.size()))
             {
-                meshPoints.append(curPoints[pointI]);
+                meshPoints.append(curPoints[pointi]);
             }
         }
     }
@@ -124,14 +121,14 @@ calcMeshData() const
     localFacesPtr_ = new List<Face>(*this);
     List<Face>& lf = *localFacesPtr_;
 
-    forAll (*this, faceI)
+    forAll(*this, facei)
     {
-        const Face& curFace = this->operator[](faceI);
-        lf[faceI].setSize(curFace.size());
+        const Face& curFace = this->operator[](facei);
+        lf[facei].setSize(curFace.size());
 
-        forAll (curFace, labelI)
+        forAll(curFace, labelI)
         {
-            lf[faceI][labelI] = markedPoints.find(curFace[labelI])();
+            lf[facei][labelI] = markedPoints.find(curFace[labelI])();
         }
     }
 
@@ -152,7 +149,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcMeshPointMap() const
@@ -182,7 +178,7 @@ calcMeshPointMap() const
     meshPointMapPtr_ = new Map<label>(2*mp.size());
     Map<label>& mpMap = *meshPointMapPtr_;
 
-    forAll (mp, i)
+    forAll(mp, i)
     {
         mpMap.insert(mp[i], i);
     }
@@ -204,7 +200,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcLocalPoints() const
@@ -235,9 +230,9 @@ calcLocalPoints() const
 
     Field<PointType>& locPts = *localPointsPtr_;
 
-    forAll (meshPts, pointI)
+    forAll(meshPts, pointi)
     {
-        locPts[pointI] = points_[meshPts[pointI]];
+        locPts[pointi] = points_[meshPts[pointi]];
     }
 
     if (debug)
@@ -257,7 +252,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcPointNormals() const
@@ -294,15 +288,15 @@ calcPointNormals() const
 
     Field<PointType>& n = *pointNormalsPtr_;
 
-    forAll (pf, pointI)
+    forAll(pf, pointi)
     {
-        PointType& curNormal = n[pointI];
+        PointType& curNormal = n[pointi];
 
-        const labelList& curFaces = pf[pointI];
+        const labelList& curFaces = pf[pointi];
 
-        forAll (curFaces, faceI)
+        forAll(curFaces, facei)
         {
-            curNormal += faceUnitNormals[curFaces[faceI]];
+            curNormal += faceUnitNormals[curFaces[facei]];
         }
 
         curNormal /= mag(curNormal) + VSMALL;
@@ -325,7 +319,56 @@ template
     class PointField,
     class PointType
 >
+void
+Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
+calcFaceCentres() const
+{
+    if (debug)
+    {
+        Pout<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
+               "calcFaceCentres() : "
+               "calculating faceCentres in PrimitivePatch"
+            << endl;
+    }
 
+    // It is considered an error to attempt to recalculate faceCentres
+    // if they have already been calculated.
+    if (faceCentresPtr_)
+    {
+        FatalErrorIn
+        (
+            "PrimitivePatch<Face, FaceList, PointField, PointType>::"
+            "calcFaceCentres()"
+        )   << "faceCentresPtr_already allocated"
+            << abort(FatalError);
+    }
+
+    faceCentresPtr_ = new Field<PointType>(this->size());
+
+    Field<PointType>& c = *faceCentresPtr_;
+
+    forAll(c, facei)
+    {
+        c[facei] = this->operator[](facei).centre(points_);
+    }
+
+    if (debug)
+    {
+        Pout<< "PrimitivePatch<Face, FaceList, PointField, PointType>::"
+               "calcFaceCentres() : "
+               "finished calculating faceCentres in PrimitivePatch"
+            << endl;
+    }
+}
+
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcFaceNormals() const
@@ -354,10 +397,10 @@ calcFaceNormals() const
 
     Field<PointType>& n = *faceNormalsPtr_;
 
-    forAll (n, faceI)
+    forAll(n, facei)
     {
-        n[faceI] = this->operator[](faceI).normal(points_);
-        n[faceI] /= mag(n[faceI]) + VSMALL;
+        n[facei] = this->operator[](facei).normal(points_);
+        n[facei] /= mag(n[facei]) + VSMALL;
     }
 
     if (debug)
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C
index 79bac6f1c0..fe65b515ec 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchMeshEdges.C
@@ -22,8 +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
-
 \*---------------------------------------------------------------------------*/
 
 #include "PrimitivePatch.H"
@@ -38,7 +36,6 @@ template
     class PointField,
     class PointType
 >
-
 Foam::labelList
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 meshEdges
@@ -116,7 +113,6 @@ template
     class PointField,
     class PointType
 >
-
 Foam::labelList
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 meshEdges
@@ -174,7 +170,6 @@ template
     class PointField,
     class PointType
 >
-
 Foam::label
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 whichEdge
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C
index 7ac3612286..e28cfc4d23 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchPointAddressing.C
@@ -40,7 +40,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcPointEdges() const
@@ -112,7 +111,6 @@ template
     class PointField,
     class PointType
 >
-
 void
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
 calcPointFaces() const
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchProjectPoints.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchProjectPoints.C
index d79caa9da4..e90bbd96ec 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchProjectPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchProjectPoints.C
@@ -42,7 +42,6 @@ template
     class PointField,
     class PointType
 >
-
 template <class ToPatch>
 Foam::List<Foam::objectHit>
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
@@ -290,7 +289,6 @@ template
     class PointField,
     class PointType
 >
-
 template <class ToPatch>
 Foam::List<Foam::objectHit>
 Foam::PrimitivePatch<Face, FaceList, PointField, PointType>::
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
index fc405b2d0e..d3f1cfe797 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C
@@ -1844,6 +1844,8 @@ void Foam::polyTopoChange::reorderCoupledFaces
     // Rotation on new faces.
     labelList rotation(faces_.size(), 0);
 
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
     // Send ordering
     forAll(boundary, patchI)
     {
@@ -1851,6 +1853,7 @@ void Foam::polyTopoChange::reorderCoupledFaces
         {
             boundary[patchI].initOrder
             (
+                pBufs,
                 primitivePatch
                 (
                     SubList<face>
@@ -1865,6 +1868,8 @@ void Foam::polyTopoChange::reorderCoupledFaces
         }
     }
 
+    pBufs.finishedSends();
+
     // Receive and calculate ordering
 
     bool anyChanged = false;
@@ -1878,6 +1883,7 @@ void Foam::polyTopoChange::reorderCoupledFaces
 
             bool changed = boundary[patchI].order
             (
+                pBufs,
                 primitivePatch
                 (
                     SubList<face>
diff --git a/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C
index e7014038ee..1b6a9e6d5a 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C
@@ -100,7 +100,7 @@ angularOscillatingDisplacementPointPatchVectorField
     angle0_(ptf.angle0_),
     amplitude_(ptf.amplitude_),
     omega_(ptf.omega_),
-    p0_(ptf.p0_)
+    p0_(ptf.p0_, mapper)
 {}
 
 
diff --git a/src/fvMotionSolver/pointPatchFields/derived/oscillatingVelocity/oscillatingVelocityPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/oscillatingVelocity/oscillatingVelocityPointPatchVectorField.C
index dbd769c557..412c9561a2 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/oscillatingVelocity/oscillatingVelocityPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/oscillatingVelocity/oscillatingVelocityPointPatchVectorField.C
@@ -91,7 +91,7 @@ oscillatingVelocityPointPatchVectorField
     fixedValuePointPatchField<vector>(ptf, p, iF, mapper),
     amplitude_(ptf.amplitude_),
     omega_(ptf.omega_),
-    p0_(ptf.p0_)
+    p0_(ptf.p0_, mapper)
 {}
 
 
diff --git a/src/fvMotionSolver/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C b/src/fvMotionSolver/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
index db0ff59bfe..8fb0408e67 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/sixDoFRigidBodyDisplacement/sixDoFRigidBodyDisplacementPointPatchVectorField.C
@@ -93,7 +93,7 @@ sixDoFRigidBodyDisplacementPointPatchVectorField
 :
     fixedValuePointPatchField<vector>(ptf, p, iF, mapper),
     motion_(ptf.motion_),
-    p0_(ptf.p0_),
+    p0_(ptf.p0_, mapper),
     rhoInf_(ptf.rhoInf_)
 {}
 
diff --git a/src/lagrangian/basic/Cloud/Cloud.C b/src/lagrangian/basic/Cloud/Cloud.C
index 5ffa5ecdc7..91e095a94a 100644
--- a/src/lagrangian/basic/Cloud/Cloud.C
+++ b/src/lagrangian/basic/Cloud/Cloud.C
@@ -105,20 +105,15 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
     const globalMeshData& pData = polyMesh_.globalData();
     const labelList& processorPatches = pData.processorPatches();
     const labelList& processorPatchIndices = pData.processorPatchIndices();
-    const labelList& processorPatchNeighbours =
-        pData.processorPatchNeighbours();
 
-    // Initialise the setpFraction moved for the particles
+    // Initialise the stepFraction moved for the particles
     forAllIter(typename Cloud<ParticleType>, *this, pIter)
     {
         pIter().stepFraction() = 0;
     }
 
-    // Assume there will be particles to transfer
-    bool transfered = true;
-
     // While there are particles to transfer
-    while (transfered)
+    while (true)
     {
         // List of lists of particles to be transfered for all the processor
         // patches
@@ -158,104 +153,93 @@ void Foam::Cloud<ParticleType>::move(TrackingData& td)
             }
         }
 
-        if (Pstream::parRun())
+        if (!Pstream::parRun())
         {
-            // List of the numbers of particles to be transfered across the
-            // processor patches
-            labelList nsTransPs(transferList.size());
+            break;
+        }
 
-            forAll(transferList, i)
+
+        // Allocate transfer buffers
+        PstreamBuffers pBufs(Pstream::nonBlocking);
+
+        // Stream into send buffers
+        forAll(transferList, i)
+        {
+            if (transferList[i].size())
             {
-                nsTransPs[i] = transferList[i].size();
-            }
-
-            // List of the numbers of particles to be transfered across the
-            // processor patches for all the processors
-            labelListList allNTrans(Pstream::nProcs());
-            allNTrans[Pstream::myProcNo()] = nsTransPs;
-            combineReduce(allNTrans, UPstream::listEq());
-
-            transfered = false;
-
-            forAll(allNTrans, i)
-            {
-                forAll(allNTrans[i], j)
-                {
-                    if (allNTrans[i][j])
-                    {
-                        transfered = true;
-                        break;
-                    }
-                }
-            }
-
-            if (!transfered)
-            {
-                break;
-            }
-
-            forAll(transferList, i)
-            {
-                if (transferList[i].size())
-                {
-                    OPstream particleStream
-                    (
-                        Pstream::blocking,
-                        refCast<const processorPolyPatch>
-                        (
-                            pMesh().boundaryMesh()[processorPatches[i]]
-                        ).neighbProcNo()
-                    );
-
-                    particleStream << transferList[i];
-                }
-            }
-
-            forAll(processorPatches, i)
-            {
-                label patchi = processorPatches[i];
-
-                const processorPolyPatch& procPatch =
+                UOPstream particleStream
+                (
                     refCast<const processorPolyPatch>
-                    (pMesh().boundaryMesh()[patchi]);
+                    (
+                        pMesh().boundaryMesh()[processorPatches[i]]
+                    ).neighbProcNo(),
+                    pBufs
+                );
 
-                label neighbProci =
-                    procPatch.neighbProcNo() - Pstream::masterNo();
+                particleStream << transferList[i];
+            }
+        }
 
-                label neighbProcPatchi = processorPatchNeighbours[patchi];
+        // Set up transfers when in non-blocking mode. Returns sizes (in bytes)
+        // to be sent/received.
+        labelListList allNTrans(Pstream::nProcs());
+        pBufs.finishedSends(allNTrans);
 
-                label nRecPs = allNTrans[neighbProci][neighbProcPatchi];
+        bool transfered = false;
 
-                if (nRecPs)
+        forAll(allNTrans, i)
+        {
+            forAll(allNTrans[i], j)
+            {
+                if (allNTrans[i][j])
                 {
-                    IPstream particleStream
-                    (
-                        Pstream::blocking,
-                        procPatch.neighbProcNo()
-                    );
-                    IDLList<ParticleType> newParticles
-                    (
-                        particleStream,
-                        typename ParticleType::iNew(*this)
-                    );
-
-                    forAllIter
-                    (
-                        typename Cloud<ParticleType>,
-                        newParticles,
-                        newpIter
-                    )
-                    {
-                        ParticleType& newp = newpIter();
-                        newp.correctAfterParallelTransfer(patchi, td);
-                        addParticle(newParticles.remove(&newp));
-                    }
+                    transfered = true;
+                    break;
                 }
             }
         }
-        else
+
+        if (!transfered)
         {
-            transfered = false;
+            break;
+        }
+
+
+        // Retrieve from receive buffers
+        forAll(processorPatches, i)
+        {
+            label patchi = processorPatches[i];
+
+            const processorPolyPatch& procPatch =
+                refCast<const processorPolyPatch>
+                (pMesh().boundaryMesh()[patchi]);
+
+            label neighbProci = procPatch.neighbProcNo();
+
+            label nRecPs = allNTrans[neighbProci][Pstream::myProcNo()];
+
+            if (nRecPs)
+            {
+                UIPstream particleStream(neighbProci, pBufs);
+
+                IDLList<ParticleType> newParticles
+                (
+                    particleStream,
+                    typename ParticleType::iNew(*this)
+                );
+
+                forAllIter
+                (
+                    typename Cloud<ParticleType>,
+                    newParticles,
+                    newpIter
+                )
+                {
+                    ParticleType& newp = newpIter();
+                    newp.correctAfterParallelTransfer(patchi, td);
+                    addParticle(newParticles.remove(&newp));
+                }
+            }
         }
     }
 }
diff --git a/src/lagrangian/dsmc/Make/files b/src/lagrangian/dsmc/Make/files
index 9f6d33d209..983de123ab 100644
--- a/src/lagrangian/dsmc/Make/files
+++ b/src/lagrangian/dsmc/Make/files
@@ -4,9 +4,6 @@ parcels/derived/dsmcParcel/dsmcParcel.C
 /* Cloud base classes */
 clouds/baseClasses/DsmcBaseCloud/DsmcBaseCloud.C
 
-/* Clouds */
-clouds/derived/dsmcCloud/dsmcCloud.C
-
 /* submodels */
 parcels/derived/dsmcParcel/defineDsmcParcel.C
 parcels/derived/dsmcParcel/makeDsmcParcelBinaryCollisionModels.C
diff --git a/src/lagrangian/dsmc/Make/options b/src/lagrangian/dsmc/Make/options
index 15874b7b55..55865dfabc 100644
--- a/src/lagrangian/dsmc/Make/options
+++ b/src/lagrangian/dsmc/Make/options
@@ -1,7 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 LIB_LIBS = \
     -llagrangian \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lmeshTools
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
index 5ac29101e5..f7aaaa724d 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.C
@@ -297,7 +297,7 @@ void Foam::DsmcCloud<ParcelType>::collisions()
     // Temporary storage for subCells
     List<DynamicList<label> > subCells(8);
 
-    scalar deltaT = cachedDeltaT();
+    scalar deltaT = mesh().time().deltaTValue();
 
     label collisionCandidates = 0;
 
@@ -473,21 +473,92 @@ void Foam::DsmcCloud<ParcelType>::collisions()
 
 
 template<class ParcelType>
-void Foam::DsmcCloud<ParcelType>::resetSurfaceDataFields()
+void Foam::DsmcCloud<ParcelType>::resetFields()
 {
-    volScalarField::GeometricBoundaryField& qBF(q_.boundaryField());
+    q_ = dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0);
 
-    forAll(qBF, p)
+    fD_ = dimensionedVector
+    (
+        "zero",
+        dimensionSet(1, -1, -2, 0, 0),
+        vector::zero
+    );
+
+    rhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
+
+    rhoM_ =  dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL);
+
+    dsmcRhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0);
+
+    linearKE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
+
+    internalE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
+
+    iDof_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
+
+    momentum_ = dimensionedVector
+    (
+        "zero",
+        dimensionSet(1, -2, -1, 0, 0),
+        vector::zero
+    );
+}
+
+
+template<class ParcelType>
+void Foam::DsmcCloud<ParcelType>::calculateFields()
+{
+    scalarField& rhoN = rhoN_.internalField();
+
+    scalarField& rhoM = rhoM_.internalField();
+
+    scalarField& dsmcRhoN = dsmcRhoN_.internalField();
+
+    scalarField& linearKE = linearKE_.internalField();
+
+    scalarField& internalE = internalE_.internalField();
+
+    scalarField& iDof = iDof_.internalField();
+
+    vectorField& momentum = momentum_.internalField();
+
+    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
     {
-        qBF[p] = 0.0;
+        const ParcelType& p = iter();
+        const label cellI = p.cell();
+
+        rhoN[cellI]++;
+
+        rhoM[cellI] += constProps(p.typeId()).mass();
+
+        dsmcRhoN[cellI]++;
+
+        linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
+
+        internalE[cellI] += p.Ei();
+
+        iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
+
+        momentum[cellI] += constProps(p.typeId()).mass()*p.U();
     }
 
-    volVectorField::GeometricBoundaryField& fDBF(fD_.boundaryField());
+    rhoN *= nParticle_/mesh().cellVolumes();
+    rhoN_.correctBoundaryConditions();
 
-    forAll(fDBF, p)
-    {
-        fDBF[p] = vector::zero;
-    }
+    rhoM *= nParticle_/mesh().cellVolumes();
+    rhoM_.correctBoundaryConditions();
+
+    linearKE *= nParticle_/mesh().cellVolumes();
+    linearKE_.correctBoundaryConditions();
+
+    internalE *= nParticle_/mesh().cellVolumes();
+    internalE_.correctBoundaryConditions();
+
+    iDof *= nParticle_/mesh().cellVolumes();
+    iDof_.correctBoundaryConditions();
+
+    momentum *= nParticle_/mesh().cellVolumes();
+    momentum_.correctBoundaryConditions();
 }
 
 
@@ -523,15 +594,14 @@ template<class ParcelType>
 Foam::DsmcCloud<ParcelType>::DsmcCloud
 (
     const word& cloudName,
-    const volScalarField& T,
-    const volVectorField& U,
+    const fvMesh& mesh,
     bool readFields
 )
 :
-    Cloud<ParcelType>(T.mesh(), cloudName, false),
+    Cloud<ParcelType>(mesh, cloudName, false),
     DsmcBaseCloud(),
     cloudName_(cloudName),
-    mesh_(T.mesh()),
+    mesh_(mesh),
     particleProperties_
     (
         IOobject
@@ -563,37 +633,142 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
     (
         IOobject
         (
-            this->name() + "q_",
+            "q",
             mesh_.time().timeName(),
             mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
         ),
-        mesh_,
-        dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
+        mesh_
     ),
     fD_
     (
         IOobject
         (
-            this->name() + "fD_",
+            "fD",
             mesh_.time().timeName(),
             mesh_,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
         ),
-        mesh_,
-        dimensionedVector
+        mesh_
+    ),
+    rhoN_
+    (
+        IOobject
         (
-            "zero",
-            dimensionSet(1, -1, -2, 0, 0),
-            vector::zero
-        )
+            "rhoN",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    rhoM_
+    (
+        IOobject
+        (
+            "rhoM",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    dsmcRhoN_
+    (
+        IOobject
+        (
+            "dsmcRhoN",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    linearKE_
+    (
+        IOobject
+        (
+            "linearKE",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    internalE_
+    (
+        IOobject
+        (
+            "internalE",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    iDof_
+    (
+        IOobject
+        (
+            "iDof",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
+    ),
+    momentum_
+    (
+        IOobject
+        (
+            "momentum",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh_
     ),
     constProps_(),
     rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
-    T_(T),
-    U_(U),
+    boundaryT_
+    (
+        volScalarField
+        (
+            IOobject
+            (
+                "boundaryT",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_
+        )
+    ),
+    boundaryU_
+    (
+        volVectorField
+        (
+            IOobject
+            (
+                "boundaryU",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            mesh_
+        )
+    ),
     binaryCollisionModel_
     (
         BinaryCollisionModel<DsmcCloud<ParcelType> >::New
@@ -641,7 +816,8 @@ template<class ParcelType>
 Foam::DsmcCloud<ParcelType>::DsmcCloud
 (
     const word& cloudName,
-    const fvMesh& mesh
+    const fvMesh& mesh,
+    const IOdictionary& dsmcInitialiseDict
 )
     :
     Cloud<ParcelType>(mesh, cloudName, false),
@@ -707,15 +883,111 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
             vector::zero
         )
     ),
+    rhoN_
+    (
+        IOobject
+        (
+            this->name() + "rhoN_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL)
+    ),
+    rhoM_
+    (
+        IOobject
+        (
+            this->name() + "rhoM_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL)
+    ),
+    dsmcRhoN_
+    (
+        IOobject
+        (
+            this->name() + "dsmcRhoN_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0)
+    ),
+    linearKE_
+    (
+        IOobject
+        (
+            this->name() + "linearKE_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0)
+    ),
+    internalE_
+    (
+        IOobject
+        (
+            this->name() + "internalE_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0)
+    ),
+    iDof_
+    (
+        IOobject
+        (
+            this->name() + "iDof_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL)
+    ),
+    momentum_
+    (
+        IOobject
+        (
+            this->name() + "momentum_",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedVector
+        (
+            "zero",
+            dimensionSet(1, -2, -1, 0, 0),
+            vector::zero
+        )
+    ),
     constProps_(),
     rndGen_(label(971501) + 1526*Pstream::myProcNo()),
-    T_
+    boundaryT_
     (
         volScalarField
         (
             IOobject
             (
-                "T",
+                "boundaryT",
                 mesh_.time().timeName(),
                 mesh_,
                 IOobject::NO_READ,
@@ -725,13 +997,13 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
             dimensionedScalar("zero",  dimensionSet(0, 0, 0, 1, 0), 0.0)
         )
     ),
-    U_
+    boundaryU_
     (
         volVectorField
         (
             IOobject
             (
-                "U",
+                "boundaryU",
                 mesh_.time().timeName(),
                 mesh_,
                 IOobject::NO_READ,
@@ -754,18 +1026,6 @@ Foam::DsmcCloud<ParcelType>::DsmcCloud
 
     buildConstProps();
 
-    IOdictionary dsmcInitialiseDict
-    (
-        IOobject
-        (
-            "dsmcInitialiseDict",
-            mesh_.time().system(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
     initialise(dsmcInitialiseDict);
 }
 
@@ -782,13 +1042,10 @@ Foam::DsmcCloud<ParcelType>::~DsmcCloud()
 template<class ParcelType>
 void Foam::DsmcCloud<ParcelType>::evolve()
 {
-    // cache the value of deltaT for this timestep
-    storeDeltaT();
-
     typename ParcelType::trackData td(*this);
 
-    // Reset the surface data collection fields
-    resetSurfaceDataFields();
+    // Reset the data collection fields
+    resetFields();
 
     if (debug)
     {
@@ -803,6 +1060,9 @@ void Foam::DsmcCloud<ParcelType>::evolve()
 
     // Calculate new velocities via stochastic collisions
     collisions();
+
+    // Calculate the volume field data
+    calculateFields();
 }
 
 
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
index c3765ba0df..5b2f382eb6 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloud.H
@@ -110,24 +110,41 @@ class DsmcCloud
         //- Force density at surface field
         volVectorField fD_;
 
+        //- number density field
+        volScalarField rhoN_;
+
+        //- Mass density field
+        volScalarField rhoM_;
+
+        //- dsmc particle density field
+        volScalarField dsmcRhoN_;
+
+        //- linear kinetic energy density field
+        volScalarField linearKE_;
+
+        //- Internal energy density field
+        volScalarField internalE_;
+
+        // Internal degree of freedom density field
+        volScalarField iDof_;
+
+        //- Momentum density field
+        volVectorField momentum_;
+
         //- Parcel constant properties - one for each type
         List<typename ParcelType::constantProperties> constProps_;
 
         //- Random number generator
         Random rndGen_;
 
-        //- In-cloud cache of deltaT, lookup in submodels and parcel is
-        //  expensive
-        scalar cachedDeltaT_;
 
+        // boundary value fields
 
-        // References to the macroscopic fields
+            //- boundary temperature
+            volScalarField boundaryT_;
 
-            //- Temperature
-            const volScalarField& T_;
-
-            //- Velocity
-            const volVectorField& U_;
+            //- boundary velocity
+            volVectorField boundaryU_;
 
 
         // References to the cloud sub-models
@@ -159,8 +176,11 @@ class DsmcCloud
         //- Calculate collisions between molecules
         void collisions();
 
-        //- Reset the surface data accumulation field values
-        void resetSurfaceDataFields();
+        //- Reset the data accumulation field values to zero
+        void resetFields();
+
+        //- Calculate the volume field data
+        void calculateFields();
 
         //- Disallow default bitwise copy construct
         DsmcCloud(const DsmcCloud&);
@@ -173,20 +193,21 @@ public:
 
     // Constructors
 
-        //- Construct given name and mesh, will read Parcels from file
+        //- Construct given name and mesh, will read Parcels and fields from
+        //  file
         DsmcCloud
         (
             const word& cloudName,
-            const volScalarField& T,
-            const volVectorField& U,
+            const fvMesh& mesh,
             bool readFields = true
         );
 
-        //- Construct given name and mesh.  Used to initialise.
+        //- Construct given name, mesh and initialisation dictionary.
         DsmcCloud
         (
             const word& cloudName,
-            const fvMesh& mesh
+            const fvMesh& mesh,
+            const IOdictionary& dsmcInitialiseDict
         );
 
 
@@ -242,35 +263,72 @@ public:
                 //- Return refernce to the random object
                 inline Random& rndGen();
 
-                //- Store (cache) the current value of deltaT
-                inline void storeDeltaT();
 
-                //- Return the cached value of deltaT
-                inline scalar cachedDeltaT() const;
+            // References to the boundary fields for surface data collection
 
+                //- Return non-const heat flux boundary field reference
+                inline volScalarField::GeometricBoundaryField& qBF();
 
-            // References to the surface data collection fields
+                //- Return non-const force density at boundary field reference
+                inline volVectorField::GeometricBoundaryField& fDBF();
 
-                //- Return heat flux at surface field
-                inline const volScalarField& q() const;
+                //- Return non-const number density boundary field reference
+                inline volScalarField::GeometricBoundaryField& rhoNBF();
 
-                //- Return non-const heat flux at surface field
-                inline volScalarField& q();
+                //- Return non-const mass density boundary field reference
+                inline volScalarField::GeometricBoundaryField& rhoMBF();
 
-                //- Return force density at surface field
-                inline const volVectorField& fD() const;
+                //- Return non-const linear kinetic energy density boundary
+                //  field reference
+                inline volScalarField::GeometricBoundaryField& linearKEBF();
 
-                //- Return non-const force density at surface field
-                inline volVectorField& fD();
+                //- Return non-const internal energy density boundary field
+                // reference
+                inline volScalarField::GeometricBoundaryField& internalEBF();
+
+                //- Return non-const internal degree of freedom density boundary
+                //  field reference
+                inline volScalarField::GeometricBoundaryField& iDofBF();
+
+                //- Return non-const momentum density boundary field reference
+                inline volVectorField::GeometricBoundaryField& momentumBF();
 
 
             // References to the macroscopic fields
 
                 //- Return macroscopic temperature
-                inline const volScalarField& T() const;
+                inline const volScalarField& boundaryT() const;
 
                 //- Return macroscopic velocity
-                inline const volVectorField& U() const;
+                inline const volVectorField& boundaryU() const;
+
+                //- Return heat flux at surface field
+                inline const volScalarField& q() const;
+
+                //- Return force density at surface field
+                inline const volVectorField& fD() const;
+
+                //- Return the real particle number density field
+                inline const volScalarField& rhoN() const;
+
+                //- Return the particle mass density field
+                inline const volScalarField& rhoM() const;
+
+                //- Return the field of number of DSMC particles
+                inline const volScalarField& dsmcRhoN() const;
+
+                //- Return the total linear kinetic energy (translational and
+                // thermal density field
+                inline const volScalarField& linearKE() const;
+
+                //- Return the internal energy density field
+                inline const volScalarField& internalE() const;
+
+                //- Return the average internal degrees of freedom  field
+                inline const volScalarField& iDof() const;
+
+                //- Return the momentum density field
+                inline const volVectorField& momentum() const;
 
 
             // Kinetic theory helper functions
@@ -385,29 +443,6 @@ public:
             void dumpParticlePositions() const;
 
 
-            // Fields
-
-                //- Return the real particle number density field
-                inline const tmp<volScalarField> rhoN() const;
-
-                //- Return the particle mass density field
-                inline const tmp<volScalarField> rhoM() const;
-
-                //- Return the field of number of DSMC particles
-                inline const tmp<volScalarField> dsmcRhoN() const;
-
-                //- Return the momentum density field
-                inline const tmp<volVectorField> momentum() const;
-
-                //- Return the total linear kinetic energy (translational and
-                // thermal density field
-                inline const tmp<volScalarField> linearKE() const;
-
-                //- Return the internal energy density field
-                inline const tmp<volScalarField> internalE() const;
-
-                //- Return the average internal degrees of freedom  field
-                inline const tmp<volScalarField> iDof() const;
 
 
         // Cloud evolution functions
diff --git a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
index 9b65d77659..beaacb2276 100644
--- a/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
+++ b/src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudI.H
@@ -126,58 +126,82 @@ inline Foam::Random& Foam::DsmcCloud<ParcelType>::rndGen()
 
 
 template<class ParcelType>
-inline void Foam::DsmcCloud<ParcelType>::storeDeltaT()
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::qBF()
 {
-    cachedDeltaT_ = mesh().time().deltaTValue();
+    return q_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline Foam::scalar Foam::DsmcCloud<ParcelType>::cachedDeltaT() const
+inline Foam::volVectorField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::fDBF()
 {
-    return cachedDeltaT_;
+    return fD_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::rhoNBF()
 {
-    return q_;
+    return rhoN_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q()
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::rhoMBF()
 {
-    return q_;
+    return rhoM_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() const
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::linearKEBF()
 {
-    return fD_;
+    return linearKE_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD()
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::internalEBF()
 {
-    return fD_;
+    return internalE_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::T() const
+inline Foam::volScalarField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::iDofBF()
 {
-    return T_;
+    return iDof_.boundaryField();
 }
 
 
 template<class ParcelType>
-inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::U() const
+inline Foam::volVectorField::GeometricBoundaryField&
+Foam::DsmcCloud<ParcelType>::momentumBF()
 {
-    return U_;
+    return momentum_.boundaryField();
+}
+
+
+template<class ParcelType>
+inline const Foam::volScalarField&
+Foam::DsmcCloud<ParcelType>::boundaryT() const
+{
+    return boundaryT_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volVectorField&
+Foam::DsmcCloud<ParcelType>::boundaryU() const
+{
+    return boundaryU_;
 }
 
 
@@ -381,265 +405,70 @@ Foam::DsmcCloud<ParcelType>::maxwellianMostProbableSpeed
 
 
 template<class ParcelType>
-inline const Foam::tmp<Foam::volScalarField>
+inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::q() const
+{
+    return q_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::fD() const
+{
+    return fD_;
+}
+
+
+template<class ParcelType>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::rhoN() const
 {
-    tmp<volScalarField> trhoN
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                this->name() + "rhoN",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
-        )
-    );
-
-    scalarField& rhoN = trhoN().internalField();
-    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
-    {
-        const ParcelType& p = iter();
-        const label cellI = p.cell();
-
-        rhoN[cellI]++;
-    }
-
-    rhoN *= nParticle_/mesh().cellVolumes();
-
-    return trhoN;
+    return rhoN_;
 }
 
 
 template<class ParcelType>
-inline const Foam::tmp<Foam::volScalarField>
-Foam::DsmcCloud<ParcelType>::rhoM() const
+inline const Foam::volScalarField& Foam::DsmcCloud<ParcelType>::rhoM() const
 {
-    tmp<volScalarField> trhoM
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                this->name() + "rhoM",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL)
-        )
-    );
-
-    scalarField& rhoM = trhoM().internalField();
-    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
-    {
-        const ParcelType& p = iter();
-        const label cellI = p.cell();
-
-        rhoM[cellI] += constProps(p.typeId()).mass();
-    }
-
-    rhoM *= nParticle_/mesh().cellVolumes();
-
-    return trhoM;
+    return rhoM_;
 }
 
 
 template<class ParcelType>
-inline const Foam::tmp<Foam::volScalarField>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::dsmcRhoN() const
 {
-    tmp<volScalarField> tdsmcRhoN
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                this->name() + "dsmcRhoN",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
-        )
-    );
-
-    scalarField& dsmcRhoN = tdsmcRhoN().internalField();
-    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
-    {
-        const ParcelType& p = iter();
-        const label cellI = p.cell();
-
-        dsmcRhoN[cellI]++;
-    }
-
-    return tdsmcRhoN;
+    return dsmcRhoN_;
 }
 
 
 template<class ParcelType>
-inline const Foam::tmp<Foam::volVectorField>
-Foam::DsmcCloud<ParcelType>::momentum() const
-{
-    tmp<volVectorField> tmomentum
-    (
-        new volVectorField
-        (
-            IOobject
-            (
-                this->name() + "momentum",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedVector
-            (
-                "zero",
-                dimensionSet(1, -2, -1, 0, 0),
-                vector::zero
-            )
-        )
-    );
-
-    vectorField& momentum = tmomentum().internalField();
-    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
-    {
-        const ParcelType& p = iter();
-        const label cellI = p.cell();
-
-        momentum[cellI] += constProps(p.typeId()).mass()*p.U();
-    }
-
-    momentum *= nParticle_/mesh().cellVolumes();
-
-    return tmomentum;
-}
-
-
-template<class ParcelType>
-inline const Foam::tmp<Foam::volScalarField>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::linearKE() const
 {
-    tmp<volScalarField> tlinearKE
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                this->name() + "linearKE",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
-        )
-    );
-
-    scalarField& linearKE = tlinearKE().internalField();
-    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
-    {
-        const ParcelType& p = iter();
-        const label cellI = p.cell();
-
-        linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
-    }
-
-    linearKE *= nParticle_/mesh().cellVolumes();
-
-    return tlinearKE;
+    return linearKE_;
 }
 
 
 template<class ParcelType>
-inline const Foam::tmp<Foam::volScalarField>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::internalE() const
 {
-    tmp<volScalarField> tinternalE
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                this->name() + "internalE",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
-        )
-    );
-
-    scalarField& internalE = tinternalE().internalField();
-    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
-    {
-        const ParcelType& p = iter();
-        const label cellI = p.cell();
-
-        internalE[cellI] += p.Ei();
-    }
-
-    internalE *= nParticle_/mesh().cellVolumes();
-
-    return tinternalE;
+    return internalE_;
 }
 
 
 template<class ParcelType>
-inline const Foam::tmp<Foam::volScalarField>
+inline const Foam::volScalarField&
 Foam::DsmcCloud<ParcelType>::iDof() const
 {
-    tmp<volScalarField> tiDof
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                this->name() + "iDof",
-                this->db().time().timeName(),
-                this->db(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE,
-                false
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
-        )
-    );
+    return iDof_;
+}
 
-    scalarField& iDof = tiDof().internalField();
 
-    forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
-    {
-        const ParcelType& p = iter();
-        const label cellI = p.cell();
-
-        iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
-    }
-
-    iDof *= nParticle_/mesh().cellVolumes();
-
-    return tiDof;
+template<class ParcelType>
+inline const Foam::volVectorField& Foam::DsmcCloud<ParcelType>::momentum() const
+{
+    return momentum_;
 }
 
 
diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
deleted file mode 100644
index 7f94f3e739..0000000000
--- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.C
+++ /dev/null
@@ -1,66 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-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
-
-\*---------------------------------------------------------------------------*/
-
-#include "dsmcCloud.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-    defineTypeNameAndDebug(dsmcCloud, 0);
-};
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::dsmcCloud::dsmcCloud
-(
-    const word& cloudName,
-    const volScalarField& T,
-    const volVectorField& U
-)
-:
-    DsmcCloud<dsmcParcel>(cloudName, T, U)
-{}
-
-
-Foam::dsmcCloud::dsmcCloud
-(
-    const word& cloudName,
-    const fvMesh& mesh
-)
-:
-    DsmcCloud<dsmcParcel>(cloudName, mesh)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::dsmcCloud::~dsmcCloud()
-{}
-
-
-// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
index 99343c2b3e..c0c17336e6 100644
--- a/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
+++ b/src/lagrangian/dsmc/clouds/derived/dsmcCloud/dsmcCloud.H
@@ -43,54 +43,8 @@ SourceFiles
 
 namespace Foam
 {
-
-/*---------------------------------------------------------------------------*\
-                          Class dsmcCloud Declaration
-\*---------------------------------------------------------------------------*/
-
-class dsmcCloud
-:
-    public DsmcCloud<dsmcParcel>
-{
-    // Private member functions
-
-        //- Disallow default bitwise copy construct
-        dsmcCloud(const dsmcCloud&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const dsmcCloud&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("dsmcCloud");
-
-
-    // Constructors
-
-        //- Construct from components
-        dsmcCloud
-        (
-            const word& cloudName,
-            const volScalarField& T,
-            const volVectorField& U
-        );
-
-        //- Construct from name and mesh, used to initialise.
-        dsmcCloud
-        (
-            const word& cloudName,
-            const fvMesh& mesh
-        );
-
-    //- Destructor
-    ~dsmcCloud();
-};
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
+    typedef DsmcCloud<dsmcParcel> dsmcCloud;
+}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
index 521559b72c..3d752aa98d 100644
--- a/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
+++ b/src/lagrangian/dsmc/parcels/Templates/DsmcParcel/DsmcParcel.C
@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "DsmcParcel.H"
+#include "meshTools.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -43,16 +44,31 @@ bool Foam::DsmcParcel<ParcelType>::move
     const polyMesh& mesh = td.cloud().pMesh();
     const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
 
-    const scalar deltaT = td.cloud().cachedDeltaT();
+    const scalar deltaT = mesh.time().deltaTValue();
     scalar tEnd = (1.0 - p.stepFraction())*deltaT;
     const scalar dtMax = tEnd;
 
+    // For reduced-D cases, the velocity used to track needs to be
+    // constrained, but the actual U_ of the parcel must not be
+    // altered or used, as it is altered by patch interactions an
+    // needs to retain its 3D value for collision purposes.
+    vector Utracking = U_;
+
     while (td.keepParticle && !td.switchProcessor && tEnd > ROOTVSMALL)
     {
+        // Apply correction to position for reduced-D cases
+        meshTools::constrainToMeshCentre(mesh, p.position());
+
+        Utracking = U_;
+
+        // Apply correction to velocity to constrain tracking for
+        // reduced-D cases
+        meshTools::constrainDirection(mesh, mesh.solutionD(), Utracking);
+
         // Set the Lagrangian time-step
         scalar dt = min(dtMax, tEnd);
 
-        dt *= p.trackToFace(p.position() + dt*U_, td);
+        dt *= p.trackToFace(p.position() + dt*Utracking, td);
 
         tEnd -= dt;
 
@@ -113,10 +129,41 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
     TrackData& td
 )
 {
+    label wppIndex = wpp.index();
+
+    label wppLocalFace = wpp.whichFace(this->face());
+
+    const scalar fA = mag(wpp.faceAreas()[wppLocalFace]);
+
+    const scalar deltaT = td.cloud().pMesh().time().deltaTValue();
+
     const constantProperties& constProps(td.cloud().constProps(typeId_));
 
     scalar m = constProps.mass();
 
+    vector nw = wpp.faceAreas()[wppLocalFace];
+    nw /= mag(nw);
+
+    scalar U_dot_nw = U_ & nw;
+
+    vector Ut = U_ - U_dot_nw*nw;
+
+    scalar invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL);
+
+    td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;
+
+    td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;
+
+    td.cloud().linearKEBF()[wppIndex][wppLocalFace] +=
+        0.5*m*(U_ & U_)*invMagUnfA;
+
+    td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;
+
+    td.cloud().iDofBF()[wppIndex][wppLocalFace] +=
+        constProps.internalDegreesOfFreedom()*invMagUnfA;
+
+    td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;
+
     // pre-interaction energy
     scalar preIE = 0.5*m*(U_ & U_) + Ei_;
 
@@ -132,27 +179,40 @@ void Foam::DsmcParcel<ParcelType>::hitWallPatch
         typeId_
     );
 
+    U_dot_nw = U_ & nw;
+
+    Ut = U_ - U_dot_nw*nw;
+
+    invMagUnfA = 1/max(mag(U_dot_nw)*fA, VSMALL);
+
+    td.cloud().rhoNBF()[wppIndex][wppLocalFace] += invMagUnfA;
+
+    td.cloud().rhoMBF()[wppIndex][wppLocalFace] += m*invMagUnfA;
+
+    td.cloud().linearKEBF()[wppIndex][wppLocalFace] +=
+        0.5*m*(U_ & U_)*invMagUnfA;
+
+    td.cloud().internalEBF()[wppIndex][wppLocalFace] += Ei_*invMagUnfA;
+
+    td.cloud().iDofBF()[wppIndex][wppLocalFace] +=
+        constProps.internalDegreesOfFreedom()*invMagUnfA;
+
+    td.cloud().momentumBF()[wppIndex][wppLocalFace] += m*Ut*invMagUnfA;
+
     // post-interaction energy
     scalar postIE = 0.5*m*(U_ & U_) + Ei_;
 
     // post-interaction momentum
     vector postIMom = m*U_;
 
-    label wppIndex = wpp.index();
-
-    label wppLocalFace = wpp.whichFace(this->face());
-
-    const scalar fA = mag(wpp.faceAreas()[wppLocalFace]);
-
-    const scalar deltaT = td.cloud().cachedDeltaT();
-
     scalar deltaQ = td.cloud().nParticle()*(preIE - postIE)/(deltaT*fA);
 
     vector deltaFD = td.cloud().nParticle()*(preIMom - postIMom)/(deltaT*fA);
 
-    td.cloud().q().boundaryField()[wppIndex][wppLocalFace] += deltaQ;
+    td.cloud().qBF()[wppIndex][wppLocalFace] += deltaQ;
+
+    td.cloud().fDBF()[wppIndex][wppLocalFace] += deltaFD;
 
-    td.cloud().fD().boundaryField()[wppIndex][wppLocalFace] += deltaFD;
 }
 
 
@@ -212,4 +272,3 @@ void Foam::DsmcParcel<ParcelType>::transformProperties
 #include "DsmcParcelIO.C"
 
 // ************************************************************************* //
-
diff --git a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C
index f15d0bdadd..b67f66ec85 100644
--- a/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C
+++ b/src/lagrangian/dsmc/parcels/derived/dsmcParcel/makeDsmcParcelWallInteractionModels.C
@@ -28,6 +28,7 @@ License
 #include "DsmcCloud.H"
 #include "MaxwellianThermal.H"
 #include "SpecularReflection.H"
+#include "MixedDiffuseSpecular.H"
 
 namespace Foam
 {
@@ -46,6 +47,12 @@ namespace Foam
         DsmcCloud,
         dsmcParcel
     );
+    makeWallInteractionModelType
+    (
+        MixedDiffuseSpecular,
+        DsmcCloud,
+        dsmcParcel
+    );
 };
 
 
diff --git a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index 03fc89d032..1df22c68d5 100644
--- a/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -129,7 +129,7 @@ void Foam::FreeStream<CloudType>::inflow()
 
     const polyMesh& mesh(cloud.mesh());
 
-    const scalar deltaT = cloud.cachedDeltaT();
+    const scalar deltaT = mesh.time().deltaTValue();
 
     Random& rndGen(cloud.rndGen());
 
@@ -139,12 +139,12 @@ void Foam::FreeStream<CloudType>::inflow()
 
     const volScalarField::GeometricBoundaryField& boundaryT
     (
-        cloud.T().boundaryField()
+        cloud.boundaryT().boundaryField()
     );
 
     const volVectorField::GeometricBoundaryField& boundaryU
     (
-        cloud.U().boundaryField()
+        cloud.boundaryU().boundaryField()
     );
 
     forAll(patches_, p)
@@ -168,7 +168,8 @@ void Foam::FreeStream<CloudType>::inflow()
             if (min(boundaryT[patchI]) < SMALL)
             {
                 FatalErrorIn ("Foam::FreeStream<CloudType>::inflow()")
-                    << "Zero boundary temperature detected, check boundaryT condition." << nl
+                    << "Zero boundary temperature detected, check boundaryT "
+                    << "condition." << nl
                     << nl << abort(FatalError);
             }
 
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index 569e9b21e6..38890400c8 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -71,10 +71,10 @@ void Foam::MaxwellianThermal<CloudType>::correct
     nw /= mag(nw);
 
     // Normal velocity magnitude
-    scalar magUn = U & nw;
+    scalar U_dot_nw = U & nw;
 
     // Wall tangential velocity (flow direction)
-    vector Ut = U - magUn*nw;
+    vector Ut = U - U_dot_nw*nw;
 
     CloudType& cloud(this->owner());
 
@@ -93,9 +93,9 @@ void Foam::MaxwellianThermal<CloudType>::correct
             U.z()*(0.8 + 0.2*rndGen.scalar01())
         );
 
-        magUn = U & nw;
+        U_dot_nw = U & nw;
 
-        Ut = U - magUn*nw;
+        Ut = U - U_dot_nw*nw;
     }
 
     // Wall tangential unit vector
@@ -104,7 +104,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
     // Other tangential unit vector
     vector tw2 = nw^tw1;
 
-    scalar T = cloud.T().boundaryField()[wppIndex][wppLocalFace];
+    scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace];
 
     scalar mass = cloud.constProps(typeId).mass();
 
@@ -118,7 +118,7 @@ void Foam::MaxwellianThermal<CloudType>::correct
           - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
         );
 
-    U += cloud.U().boundaryField()[wppIndex][wppLocalFace];
+    U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
 
     Ei = cloud.equipartitionInternalEnergy(T, iDof);
 }
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
new file mode 100644
index 0000000000..9f2f5fdc3b
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
@@ -0,0 +1,140 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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
+
+\*---------------------------------------------------------------------------*/
+
+#include "MixedDiffuseSpecular.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::MixedDiffuseSpecular<CloudType>::MixedDiffuseSpecular
+(
+    const dictionary& dict,
+    CloudType& cloud
+)
+:
+    WallInteractionModel<CloudType>(dict, cloud, typeName),
+    diffuseFraction_(readScalar(this->coeffDict().lookup("diffuseFraction")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::MixedDiffuseSpecular<CloudType>::~MixedDiffuseSpecular()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template <class CloudType>
+void Foam::MixedDiffuseSpecular<CloudType>::correct
+(
+    const wallPolyPatch& wpp,
+    const label faceId,
+    vector& U,
+    scalar& Ei,
+    label typeId
+)
+{
+    label wppIndex = wpp.index();
+
+    label wppLocalFace = wpp.whichFace(faceId);
+
+    vector nw = wpp.faceAreas()[wppLocalFace];
+
+    // Normal unit vector
+    nw /= mag(nw);
+
+    // Normal velocity magnitude
+    scalar U_dot_nw = U & nw;
+
+    CloudType& cloud(this->owner());
+
+    Random& rndGen(cloud.rndGen());
+
+    if (diffuseFraction_ > rndGen.scalar01())
+    {
+        // Diffuse reflection
+
+        // Wall tangential velocity (flow direction)
+        vector Ut = U - U_dot_nw*nw;
+
+        while (mag(Ut) < SMALL)
+        {
+            // If the incident velocity is parallel to the face normal, no
+            // tangential direction can be chosen.  Add a perturbation to the
+            // incoming velocity and recalculate.
+
+            U = vector
+            (
+                U.x()*(0.8 + 0.2*rndGen.scalar01()),
+                U.y()*(0.8 + 0.2*rndGen.scalar01()),
+                U.z()*(0.8 + 0.2*rndGen.scalar01())
+            );
+
+            U_dot_nw = U & nw;
+
+            Ut = U - U_dot_nw*nw;
+        }
+
+        // Wall tangential unit vector
+        vector tw1 = Ut/mag(Ut);
+
+        // Other tangential unit vector
+        vector tw2 = nw^tw1;
+
+        scalar T = cloud.boundaryT().boundaryField()[wppIndex][wppLocalFace];
+
+        scalar mass = cloud.constProps(typeId).mass();
+
+        scalar iDof = cloud.constProps(typeId).internalDegreesOfFreedom();
+
+        U =
+            sqrt(physicoChemical::k.value()*T/mass)
+           *(
+                rndGen.GaussNormal()*tw1
+              + rndGen.GaussNormal()*tw2
+              - sqrt(-2.0*log(max(1 - rndGen.scalar01(), VSMALL)))*nw
+            );
+
+        U += cloud.boundaryU().boundaryField()[wppIndex][wppLocalFace];
+
+        Ei = cloud.equipartitionInternalEnergy(T, iDof);
+    }
+    else
+    {
+        // Specular reflection
+
+        if (U_dot_nw > 0.0)
+        {
+            U -= 2.0*U_dot_nw*nw;
+        }
+    }
+
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H
new file mode 100644
index 0000000000..b84f43d8f6
--- /dev/null
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-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
+
+Class
+    Foam::MixedDiffuseSpecular
+
+Description
+    Wall interaction setting microscopic velocity to a random one drawn from a
+    Maxwellian distribution corresponding to a specified temperature
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MixedDiffuseSpecular_H
+#define MixedDiffuseSpecular_H
+
+#include "WallInteractionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                    Class MixedDiffuseSpecular Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class MixedDiffuseSpecular
+:
+    public WallInteractionModel<CloudType>
+{
+    // Private data
+
+        //- Fraction of wall interactions that are diffuse
+        scalar diffuseFraction_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("MixedDiffuseSpecular");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        MixedDiffuseSpecular
+        (
+            const dictionary& dict,
+            CloudType& cloud
+        );
+
+
+    // Destructor
+    virtual ~MixedDiffuseSpecular();
+
+
+    // Member Functions
+
+        //- Apply wall correction
+        virtual void correct
+        (
+            const wallPolyPatch& wpp,
+            const label faceId,
+            vector& U,
+            scalar& Ei,
+            label typeId
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "MixedDiffuseSpecular.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
index 1fb5ab3208..bed84b8c68 100644
--- a/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
+++ b/src/lagrangian/dsmc/submodels/WallInteractionModel/SpecularReflection/SpecularReflection.C
@@ -63,11 +63,11 @@ void Foam::SpecularReflection<CloudType>::correct
     vector nw = wpp.faceAreas()[wpp.whichFace(faceId)];
     nw /= mag(nw);
 
-    scalar magUn = U & nw;
+    scalar U_dot_nw = U & nw;
 
-    if (magUn > 0.0)
+    if (U_dot_nw > 0.0)
     {
-        U -= 2.0*magUn*nw;
+        U -= 2.0*U_dot_nw*nw;
     }
 }
 
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C
index 439ed64bdf..7cc9e70c0e 100644
--- a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.C
@@ -25,27 +25,49 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "distribution.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+    defineTypeNameAndDebug(distribution, 0);
+}
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+void Foam::distribution::write
+(
+    const fileName& file,
+    const List<Pair<scalar> >& pairs
+)
+{
+    OFstream os(file);
+
+    forAll(pairs, i)
+    {
+        os  << pairs[i].first() << ' ' << pairs[i].second() << nl;
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-distribution::distribution()
+Foam::distribution::distribution()
 :
     Map<label>(),
     binWidth_(1)
 {}
 
 
-distribution::distribution(const scalar binWidth)
+Foam::distribution::distribution(const scalar binWidth)
 :
     Map<label>(),
     binWidth_(binWidth)
 {}
 
 
-distribution::distribution(const distribution& d)
+Foam::distribution::distribution(const distribution& d)
 :
     Map<label>(static_cast< Map<label> >(d)),
     binWidth_(d.binWidth())
@@ -54,13 +76,13 @@ distribution::distribution(const distribution& d)
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-distribution::~distribution()
+Foam::distribution::~distribution()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-label distribution::totalEntries() const
+Foam::label Foam::distribution::totalEntries() const
 {
     label sumOfEntries = 0;
 
@@ -88,7 +110,7 @@ label distribution::totalEntries() const
 }
 
 
-scalar distribution::approxTotalEntries() const
+Foam::scalar Foam::distribution::approxTotalEntries() const
 {
     scalar sumOfEntries = 0;
 
@@ -101,7 +123,7 @@ scalar distribution::approxTotalEntries() const
 }
 
 
-scalar distribution::mean() const
+Foam::scalar Foam::distribution::mean() const
 {
     scalar runningSum = 0;
 
@@ -124,7 +146,7 @@ scalar distribution::mean() const
 }
 
 
-scalar distribution::median()
+Foam::scalar Foam::distribution::median()
 {
     // From:
     // http://mathworld.wolfram.com/StatisticalMedian.html
@@ -188,7 +210,7 @@ scalar distribution::median()
 }
 
 
-void distribution::add(const scalar valueToAdd)
+void Foam::distribution::add(const scalar valueToAdd)
 {
     iterator iter(this->begin());
 
@@ -218,13 +240,13 @@ void distribution::add(const scalar valueToAdd)
 }
 
 
-void distribution::add(const label valueToAdd)
+void Foam::distribution::add(const label valueToAdd)
 {
     add(scalar(valueToAdd));
 }
 
 
-void distribution::insertMissingKeys()
+void Foam::distribution::insertMissingKeys()
 {
     iterator iter(this->begin());
 
@@ -247,7 +269,7 @@ void distribution::insertMissingKeys()
 }
 
 
-List< Pair<scalar> > distribution::normalised()
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalised()
 {
     scalar totEnt = approxTotalEntries();
 
@@ -268,17 +290,25 @@ List< Pair<scalar> > distribution::normalised()
         normDist[k].second() = scalar((*this)[key])/totEnt/binWidth_;
     }
 
+    if (debug)
+    {
+        Info<< "totEnt: " << totEnt << endl;
+    }
+
     return normDist;
 }
 
 
-List< Pair<scalar> > distribution::normalisedMinusMean()
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalisedMinusMean()
 {
     return normalisedShifted(mean());
 }
 
 
-List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::normalisedShifted
+(
+    scalar shiftValue
+)
 {
     List<Pair<scalar> > oldDist(normalised());
 
@@ -301,20 +331,23 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
 
     label newKey = lowestNewKey;
 
-//     Info << shiftValue
-//         << nl << lowestOldBin
-//         << nl << lowestNewKey
-//         << nl << interpolationStartDirection
-//         << endl;
+    if (debug)
+    {
+        Info<< shiftValue
+            << nl << lowestOldBin
+            << nl << lowestNewKey
+            << nl << interpolationStartDirection
+            << endl;
 
-//     scalar checkNormalisation = 0;
+        scalar checkNormalisation = 0;
 
-//     forAll (oldDist, oD)
-//     {
-//         checkNormalisation += oldDist[oD].second()*binWidth_;
-//     }
+        forAll (oldDist, oD)
+        {
+            checkNormalisation += oldDist[oD].second()*binWidth_;
+        }
 
-//     Info << "Initial normalisation = " << checkNormalisation << endl;
+        Info<< "Initial normalisation = " << checkNormalisation << endl;
+    }
 
     forAll(oldDist,u)
     {
@@ -368,20 +401,23 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
         newKey++;
     }
 
-//     checkNormalisation = 0;
+    if (debug)
+    {
+        scalar checkNormalisation = 0;
 
-//     forAll (newDist, nD)
-//     {
-//         checkNormalisation += newDist[nD].second()*binWidth_;
-//     }
+        forAll (newDist, nD)
+        {
+            checkNormalisation += newDist[nD].second()*binWidth_;
+        }
 
-//     Info << "Shifted normalisation = " << checkNormalisation << endl;
+        Info<< "Shifted normalisation = " << checkNormalisation << endl;
+    }
 
     return newDist;
 }
 
 
-List<Pair<scalar> > distribution::raw()
+Foam::List<Foam::Pair<Foam::scalar> > Foam::distribution::raw()
 {
     insertMissingKeys();
 
@@ -406,7 +442,7 @@ List<Pair<scalar> > distribution::raw()
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-void distribution::operator=(const distribution& rhs)
+void Foam::distribution::operator=(const distribution& rhs)
 {
     // Check for assignment to self
     if (this == &rhs)
@@ -424,7 +460,7 @@ void distribution::operator=(const distribution& rhs)
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-Ostream& operator<<(Ostream& os, const distribution& d)
+Foam::Ostream& Foam::operator<<(Ostream& os, const distribution& d)
 {
     os  << d.binWidth_
         << static_cast<const Map<label>&>(d);
@@ -440,8 +476,4 @@ Ostream& operator<<(Ostream& os, const distribution& d)
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H
index da1ce973bc..d86ae83934 100644
--- a/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H
+++ b/src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.H
@@ -26,6 +26,8 @@ Class
     Foam::distribution
 
 Description
+    Accumulating histogram of values.  Specified bin resolution
+    automatic generation of bins.
 
 SourceFiles
     distributionI.H
@@ -59,6 +61,21 @@ class distribution
 
 public:
 
+    //- Runtime type information
+
+        TypeName("distribution");
+
+    // Static functions
+
+        //- write to file
+
+            static void write
+            (
+                const fileName& file,
+                const List<Pair<scalar> >& pairs
+            );
+
+
     // Constructors
 
         //- Construct null
@@ -73,7 +90,7 @@ public:
 
     // Destructor
 
-        ~distribution();
+        virtual ~distribution();
 
 
     // Member Functions
@@ -97,7 +114,7 @@ public:
 
         List<Pair<scalar> > normalisedMinusMean();
 
-        List<Pair<scalar> > normalisedShifted(const scalar shiftValue);
+        List<Pair<scalar> > normalisedShifted(scalar shiftValue);
 
         List<Pair<scalar> > raw();
 
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 88f806b5dd..307cb743f8 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -2411,6 +2411,9 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
     }
     syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
 
+    // Get per face whether is it master (of a coupled set of faces)
+    PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
+
     // Set owner as no-flip
     forAll(patches, patchI)
     {
@@ -2428,14 +2431,27 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
 
                 bool flip;
-                if (ownZone == max(ownZone, neiZone))
+
+                label maxZone = max(ownZone, neiZone);
+
+                if (maxZone == -1)
                 {
                     flip = false;
                 }
-                else
+                else if (ownZone == neiZone)
+                {
+                    // Can only happen for coupled boundaries. Keep master
+                    // face unflipped.
+                    flip = !isMasterFace[faceI];
+                }
+                else if (neiZone == maxZone)
                 {
                     flip = true;
                 }
+                else
+                {
+                    flip = false;
+                }
 
                 meshMod.setAction
                 (
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
index e16b8204f7..0da6068681 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.C
@@ -145,44 +145,49 @@ Foam::directMappedPolyPatch::~directMappedPolyPatch()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 //- Initialise the calculation of the patch geometry
-void Foam::directMappedPolyPatch::initGeometry()
+void Foam::directMappedPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
-    polyPatch::initGeometry();
-    directMappedPatchBase::clearOut();
+    polyPatch::initGeometry(pBufs);
 }
 
 //- Calculate the patch geometry
-void Foam::directMappedPolyPatch::calcGeometry()
+void Foam::directMappedPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 {
-    polyPatch::calcGeometry();
+    polyPatch::calcGeometry(pBufs);
     directMappedPatchBase::clearOut();
 }
 
 //- Initialise the patches for moving points
-void Foam::directMappedPolyPatch::initMovePoints(const pointField& p)
+void Foam::directMappedPolyPatch::initMovePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
 {
-    polyPatch::initMovePoints(p);
-    directMappedPatchBase::clearOut();
+    polyPatch::initMovePoints(pBufs, p);
 }
 
 //- Correct patches after moving points
-void Foam::directMappedPolyPatch::movePoints(const pointField& p)
+void Foam::directMappedPolyPatch::movePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
 {
-    polyPatch::movePoints(p);
+    polyPatch::movePoints(pBufs, p);
     directMappedPatchBase::clearOut();
 }
 
 //- Initialise the update of the patch topology
-void Foam::directMappedPolyPatch::initUpdateMesh()
+void Foam::directMappedPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    polyPatch::initUpdateMesh();
-    directMappedPatchBase::clearOut();
+    polyPatch::initUpdateMesh(pBufs);
 }
 
 //- Update of the patch topology
-void Foam::directMappedPolyPatch::updateMesh()
+void Foam::directMappedPolyPatch::updateMesh(PstreamBuffers& pBufs)
 {
-    polyPatch::updateMesh();
+    polyPatch::updateMesh(pBufs);
     directMappedPatchBase::clearOut();
 }
 
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
index 8cb5907c95..654c3470a4 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedPolyPatch.H
@@ -65,22 +65,22 @@ class directMappedPolyPatch
 protected:
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry();
+        virtual void initGeometry(PstreamBuffers&);
 
         //- Calculate the patch geometry
-        virtual void calcGeometry();
+        virtual void calcGeometry(PstreamBuffers&);
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&);
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&);
+        virtual void movePoints(PstreamBuffers&, const pointField&);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh();
+        virtual void initUpdateMesh(PstreamBuffers&);
 
         //- Update of the patch topology
-        virtual void updateMesh();
+        virtual void updateMesh(PstreamBuffers&);
 
 
 public:
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C b/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C
index 6d6038430d..c4295a31a4 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.C
@@ -150,44 +150,49 @@ Foam::directMappedWallPolyPatch::~directMappedWallPolyPatch()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 //- Initialise the calculation of the patch geometry
-void Foam::directMappedWallPolyPatch::initGeometry()
+void Foam::directMappedWallPolyPatch::initGeometry(PstreamBuffers& pBufs)
 {
-    wallPolyPatch::initGeometry();
-    directMappedPatchBase::clearOut();
+    wallPolyPatch::initGeometry(pBufs);
 }
 
 //- Calculate the patch geometry
-void Foam::directMappedWallPolyPatch::calcGeometry()
+void Foam::directMappedWallPolyPatch::calcGeometry(PstreamBuffers& pBufs)
 {
-    wallPolyPatch::calcGeometry();
+    wallPolyPatch::calcGeometry(pBufs);
     directMappedPatchBase::clearOut();
 }
 
 //- Initialise the patches for moving points
-void Foam::directMappedWallPolyPatch::initMovePoints(const pointField& p)
+void Foam::directMappedWallPolyPatch::initMovePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
 {
-    wallPolyPatch::initMovePoints(p);
-    directMappedPatchBase::clearOut();
+    wallPolyPatch::initMovePoints(pBufs, p);
 }
 
 //- Correct patches after moving points
-void Foam::directMappedWallPolyPatch::movePoints(const pointField& p)
+void Foam::directMappedWallPolyPatch::movePoints
+(
+    PstreamBuffers& pBufs,
+    const pointField& p
+)
 {
-    wallPolyPatch::movePoints(p);
+    wallPolyPatch::movePoints(pBufs, p);
     directMappedPatchBase::clearOut();
 }
 
 //- Initialise the update of the patch topology
-void Foam::directMappedWallPolyPatch::initUpdateMesh()
+void Foam::directMappedWallPolyPatch::initUpdateMesh(PstreamBuffers& pBufs)
 {
-    wallPolyPatch::initUpdateMesh();
-    directMappedPatchBase::clearOut();
+    wallPolyPatch::initUpdateMesh(pBufs);
 }
 
 //- Update of the patch topology
-void Foam::directMappedWallPolyPatch::updateMesh()
+void Foam::directMappedWallPolyPatch::updateMesh(PstreamBuffers& pBufs)
 {
-    wallPolyPatch::updateMesh();
+    wallPolyPatch::updateMesh(pBufs);
     directMappedPatchBase::clearOut();
 }
 
diff --git a/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H b/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H
index 0a2f5bdfd3..b0657164da 100644
--- a/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H
+++ b/src/meshTools/directMapped/directMappedPolyPatch/directMappedWallPolyPatch.H
@@ -65,22 +65,22 @@ class directMappedWallPolyPatch
 protected:
 
         //- Initialise the calculation of the patch geometry
-        virtual void initGeometry();
+        virtual void initGeometry(PstreamBuffers&);
 
         //- Calculate the patch geometry
-        virtual void calcGeometry();
+        virtual void calcGeometry(PstreamBuffers&);
 
         //- Initialise the patches for moving points
-        virtual void initMovePoints(const pointField&);
+        virtual void initMovePoints(PstreamBuffers&, const pointField&);
 
         //- Correct patches after moving points
-        virtual void movePoints(const pointField&);
+        virtual void movePoints(PstreamBuffers&, const pointField&);
 
         //- Initialise the update of the patch topology
-        virtual void initUpdateMesh();
+        virtual void initUpdateMesh(PstreamBuffers&);
 
         //- Update of the patch topology
-        virtual void updateMesh();
+        virtual void updateMesh(PstreamBuffers&);
 
 
 public:
diff --git a/src/postProcessing/functionObjects/forces/forces/forces.H b/src/postProcessing/functionObjects/forces/forces/forces.H
index 42bb67c23b..e3f7cf03fb 100644
--- a/src/postProcessing/functionObjects/forces/forces/forces.H
+++ b/src/postProcessing/functionObjects/forces/forces/forces.H
@@ -160,6 +160,9 @@ protected:
         //- If the forces file has not been created create it
         void makeFile();
 
+        //- Output file header information
+        virtual void writeFileHeader();
+
         //- Return the effective viscous stress (laminar + turbulent).
         tmp<volSymmTensorField> devRhoReff() const;
 
@@ -176,9 +179,6 @@ protected:
         //- Disallow default bitwise assignment
         void operator=(const forces&);
 
-        //- Output file header information
-        virtual void writeFileHeader();
-
 
 public:
 
diff --git a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
index ca2f0b3467..7d72ddd265 100644
--- a/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
+++ b/src/postProcessing/functionObjects/utilities/dsmcFields/dsmcFields.C
@@ -115,6 +115,7 @@ void Foam::dsmcFields::write()
         word linearKEMeanName = "linearKEMean";
         word internalEMeanName = "internalEMean";
         word iDofMeanName = "iDofMean";
+        word fDMeanName = "fDMean";
 
         const volScalarField& rhoNMean = obr_.lookupObject<volScalarField>
         (
@@ -146,6 +147,11 @@ void Foam::dsmcFields::write()
             iDofMeanName
         );
 
+        volVectorField fDMean =  obr_.lookupObject<volVectorField>
+        (
+            fDMeanName
+        );
+
         if (min(mag(rhoNMean)).value() > VSMALL)
         {
             Info<< "Calculating dsmcFields." << endl;
@@ -173,8 +179,9 @@ void Foam::dsmcFields::write()
                     obr_,
                     IOobject::NO_READ
                 ),
+
                 2.0/(3.0*physicoChemical::k.value()*rhoNMean)
-                *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
+               *(linearKEMean - 0.5*rhoMMean*(UMean & UMean))
             );
 
             Info<< "    Calculating internalT field." << endl;
@@ -187,7 +194,7 @@ void Foam::dsmcFields::write()
                     obr_,
                     IOobject::NO_READ
                 ),
-                2.0/(physicoChemical::k.value()*iDofMean)*internalEMean
+                (2.0/physicoChemical::k.value())*(internalEMean/iDofMean)
             );
 
             Info<< "    Calculating overallT field." << endl;
@@ -201,9 +208,36 @@ void Foam::dsmcFields::write()
                     IOobject::NO_READ
                 ),
                 2.0/(physicoChemical::k.value()*(3.0*rhoNMean + iDofMean))
-                *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
+               *(linearKEMean - 0.5*rhoMMean*(UMean & UMean) + internalEMean)
             );
 
+            Info<< "    Calculating pressure field." << endl;
+            volScalarField p
+            (
+                IOobject
+                (
+                    "p",
+                    obr_.time().timeName(),
+                    obr_,
+                    IOobject::NO_READ
+                ),
+                physicoChemical::k.value()*rhoNMean*translationalT
+            );
+
+            const fvMesh& mesh = fDMean.mesh();
+
+            forAll(mesh.boundaryMesh(), i)
+            {
+                const polyPatch& patch = mesh.boundaryMesh()[i];
+
+                if (isA<wallPolyPatch>(patch))
+                {
+                    p.boundaryField()[i] =
+                        fDMean.boundaryField()[i]
+                      & (patch.faceAreas()/mag(patch.faceAreas()));
+                }
+            }
+
             Info<< "    mag(UMean) max/min : "
                 << max(mag(UMean)).value() << " "
                 << min(mag(UMean)).value() << endl;
@@ -220,6 +254,10 @@ void Foam::dsmcFields::write()
                 << max(overallT).value() << " "
                 << min(overallT).value() << endl;
 
+            Info<< "    p max/min : "
+                << max(p).value() << " "
+                << min(p).value() << endl;
+
             UMean.write();
 
             translationalT.write();
@@ -228,6 +266,8 @@ void Foam::dsmcFields::write()
 
             overallT.write();
 
+            p.write();
+
             Info<< "dsmcFields written." << nl << endl;
         }
         else
diff --git a/src/sampling/cuttingPlane/cuttingPlane.C b/src/sampling/cuttingPlane/cuttingPlane.C
index 09585ec5e7..80a4e72bd5 100644
--- a/src/sampling/cuttingPlane/cuttingPlane.C
+++ b/src/sampling/cuttingPlane/cuttingPlane.C
@@ -258,6 +258,7 @@ bool Foam::cuttingPlane::walkCell
 void Foam::cuttingPlane::walkCellCuts
 (
     const primitiveMesh& mesh,
+    const bool triangulate,
     const UList<label>& edgePoint
 )
 {
@@ -293,7 +294,7 @@ void Foam::cuttingPlane::walkCellCuts
         // Check for the unexpected ...
         if (startEdgeI == -1)
         {
-            FatalErrorIn("Foam::cuttingPlane::walkCellCuts")
+            FatalErrorIn("Foam::cuttingPlane::walkCellCuts(..)")
                 << "Cannot find cut edge for cut cell " << cellI
                 << abort(FatalError);
         }
@@ -318,10 +319,18 @@ void Foam::cuttingPlane::walkCellCuts
                 f = f.reverseFace();
             }
 
-            // the cut faces are usually quite ugly, so always triangulate
-            label nTri = f.triangles(cutPoints, dynCutFaces);
-            while (nTri--)
+            // the cut faces are usually quite ugly, so optionally triangulate
+            if (triangulate)
             {
+                label nTri = f.triangles(cutPoints, dynCutFaces);
+                while (nTri--)
+                {
+                    dynCutCells.append(cellI);
+                }
+            }
+            else
+            {
+                dynCutFaces.append(f);
                 dynCutCells.append(cellI);
             }
         }
@@ -346,12 +355,13 @@ Foam::cuttingPlane::cuttingPlane
 (
     const plane& pln,
     const primitiveMesh& mesh,
+    const bool triangulate,
     const UList<label>& cellIdLabels
 )
 :
     plane(pln)
 {
-    reCut(mesh, cellIdLabels);
+    reCut(mesh, triangulate, cellIdLabels);
 }
 
 
@@ -362,6 +372,7 @@ Foam::cuttingPlane::cuttingPlane
 void Foam::cuttingPlane::reCut
 (
     const primitiveMesh& mesh,
+    const bool triangulate,
     const UList<label>& cellIdLabels
 )
 {
@@ -379,7 +390,7 @@ void Foam::cuttingPlane::reCut
     intersectEdges(mesh, dotProducts, edgePoint);
 
     // Do topological walk around cell to find closed loop.
-    walkCellCuts(mesh, edgePoint);
+    walkCellCuts(mesh, triangulate, edgePoint);
 }
 
 
diff --git a/src/sampling/cuttingPlane/cuttingPlane.H b/src/sampling/cuttingPlane/cuttingPlane.H
index 6804a9ab18..aba51a78ed 100644
--- a/src/sampling/cuttingPlane/cuttingPlane.H
+++ b/src/sampling/cuttingPlane/cuttingPlane.H
@@ -105,6 +105,7 @@ class cuttingPlane
         void walkCellCuts
         (
             const primitiveMesh& mesh,
+            const bool triangulate,
             const UList<label>& edgePoint
         );
 
@@ -122,6 +123,7 @@ protected:
         void reCut
         (
             const primitiveMesh&,
+            const bool triangulate,
             const UList<label>& cellIdLabels = UList<label>::null()
         );
 
@@ -138,6 +140,7 @@ public:
         (
             const plane&,
             const primitiveMesh&,
+            const bool triangulate,
             const UList<label>& cellIdLabels = UList<label>::null()
         );
 
diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
index c5ce93f1d0..1f03cdcf36 100644
--- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
+++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C
@@ -150,11 +150,11 @@ bool Foam::sampledPlane::update()
 
     if (zoneId < 0)
     {
-        reCut(mesh());
+        reCut(mesh(), true);    // always triangulate. Note:Make option?
     }
     else
     {
-        reCut(mesh(), mesh().cellZones()[zoneId]);
+        reCut(mesh(), true, mesh().cellZones()[zoneId]);
     }
 
     if (debug)
diff --git a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H
index 717f6b1e7c..c5e57745d7 100644
--- a/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H
+++ b/src/turbulenceModels/compressible/RAS/derivedFvPatchFields/turbulentTemperatureCoupledBaffle/turbulentTemperatureCoupledBaffleFvPatchScalarField.H
@@ -97,7 +97,7 @@ class turbulentTemperatureCoupledBaffleFvPatchScalarField
 public:
 
     //- Runtime type information
-    TypeName("turbulentTemperatureCoupledBaffle");
+    TypeName("compressible::turbulentTemperatureCoupledBaffle");
 
 
     // Constructors
diff --git a/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.C b/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.C
index 53a0773457..60bbb41ef1 100644
--- a/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.C
+++ b/src/turbulenceModels/incompressible/LES/Smagorinsky/Smagorinsky.C
@@ -73,6 +73,8 @@ Smagorinsky::Smagorinsky
     )
 {
     updateSubGridScaleFields(fvc::grad(U));
+
+    printCoeffs();
 }
 
 
diff --git a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/iDof b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/iDof
index 4f081fa179..1e13af9370 100644
--- a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/iDof
+++ b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/iDof
@@ -32,7 +32,8 @@ boundaryField
 
     walls
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 
 }
diff --git a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/internalE b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/internalE
index 8bf198ceae..976ec30e74 100644
--- a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/internalE
+++ b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/internalE
@@ -32,7 +32,8 @@ boundaryField
 
     walls
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 }
 
diff --git a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/linearKE b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/linearKE
index f05be3b1c8..b26e99313c 100644
--- a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/linearKE
+++ b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/linearKE
@@ -32,7 +32,8 @@ boundaryField
 
     walls
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 }
 
diff --git a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/momentum b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/momentum
index 244901de0c..05d234bd2a 100644
--- a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/momentum
+++ b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/momentum
@@ -32,7 +32,8 @@ boundaryField
 
     walls
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform (0 0 0);
     }
 
 }
diff --git a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoM b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoM
index fa4dec993b..7ec6f3368f 100644
--- a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoM
+++ b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoM
@@ -32,7 +32,8 @@ boundaryField
 
     walls
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 }
 
diff --git a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoN b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoN
index cd16296be9..f847c22cd0 100644
--- a/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoN
+++ b/tutorials/discreteMethods/dsmcFoam/supersonicCorner/0/rhoN
@@ -32,7 +32,8 @@ boundaryField
 
     walls
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 }
 
diff --git a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/iDof b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/iDof
index 3e8f63f7c3..b31b782ef0 100644
--- a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/iDof
+++ b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/iDof
@@ -27,7 +27,8 @@ boundaryField
 
     obstacle
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 
     periodic
diff --git a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/internalE b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/internalE
index 4cb2411197..61e30470b9 100644
--- a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/internalE
+++ b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/internalE
@@ -27,7 +27,8 @@ boundaryField
 
     obstacle
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 
     periodic
diff --git a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/linearKE b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/linearKE
index 72f2c152c5..b3f5a50d70 100644
--- a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/linearKE
+++ b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/linearKE
@@ -27,7 +27,8 @@ boundaryField
 
     obstacle
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 
     periodic
diff --git a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/momentum b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/momentum
index f0e6facaef..424e0310cc 100644
--- a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/momentum
+++ b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/momentum
@@ -27,7 +27,8 @@ boundaryField
 
     obstacle
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform (0 0 0);
     }
 
     periodic
diff --git a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoM b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoM
index 9a84df6fa0..224aae627a 100644
--- a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoM
+++ b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoM
@@ -27,7 +27,8 @@ boundaryField
 
     obstacle
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 
     periodic
diff --git a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoN b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoN
index 4f948784ae..275417c995 100644
--- a/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoN
+++ b/tutorials/discreteMethods/dsmcFoam/wedge15Ma5/0/rhoN
@@ -27,7 +27,8 @@ boundaryField
 
     obstacle
     {
-        type            zeroGradient;
+        type            calculated;
+        value           uniform 0;
     }
 
     periodic
diff --git a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/decomposeParDict b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/decomposeParDict
index 176b1467df..6e02e44e34 100644
--- a/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/decomposeParDict
+++ b/tutorials/lagrangian/rhoPisoTwinParcelFoam/simplifiedSiwek/system/decomposeParDict
@@ -34,7 +34,7 @@ hierarchicalCoeffs
 
 metisCoeffs
 {
-    processorWeights ( 1 1 1 1 );
+    //processorWeights ( 1 1 1 1 );
 }
 
 manualCoeffs