diff --git a/applications/utilities/mesh/conversion/gambitToFoam/gambitToFoam.L b/applications/utilities/mesh/conversion/gambitToFoam/gambitToFoam.L
index 7706ccb272..dc4531fef6 100644
--- a/applications/utilities/mesh/conversion/gambitToFoam/gambitToFoam.L
+++ b/applications/utilities/mesh/conversion/gambitToFoam/gambitToFoam.L
@@ -433,7 +433,7 @@ mtype                 {space}"MTYPE:"{space}
         }
         else
         {
-            curGroupID = readLabel(groupStream);;
+            curGroupID = readLabel(groupStream);
         }
 
         BEGIN(cellStreams);
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/Make/options b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
index 2981cf30be..dc78db3436 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/Make/options
+++ b/applications/utilities/mesh/manipulation/renumberMesh/Make/options
@@ -11,4 +11,4 @@ EXE_LIBS = \
     -lfiniteVolume \
     -lgenericPatchFields \
     -lrenumberMethods \
-    -ldecompositionMethods
+    -ldecompositionMethods -L$(FOAM_LIBBIN)/dummy -lmetisDecomp -lscotchDecomp
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
index 7e090fa30e..02e5a06567 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C
@@ -29,7 +29,7 @@ Description
     renumbering all fields from all the time directories.
 
     By default uses bandCompression (CuthillMcKee) but will
-    read system/renumberMeshDict if present and use the method from there.
+    read system/renumberMeshDict if -dict option is present
 
 \*---------------------------------------------------------------------------*/
 
@@ -45,6 +45,7 @@ Description
 #include "renumberMethod.H"
 #include "zeroGradientFvPatchFields.H"
 #include "CuthillMcKeeRenumber.H"
+#include "fvMeshSubset.H"
 
 using namespace Foam;
 
@@ -455,6 +456,69 @@ autoPtr<mapPolyMesh> reorderMesh
 }
 
 
+// Return new to old cell numbering
+labelList regionRenumber
+(
+    const renumberMethod& method,
+    const fvMesh& mesh,
+    const labelList& cellToRegion
+)
+{
+    Info<< "Determining cell order:" << endl;
+
+    labelList cellOrder(cellToRegion.size());
+
+    label nRegions = max(cellToRegion)+1;
+
+    labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
+
+    label cellI = 0;
+
+    forAll(regionToCells, regionI)
+    {
+        Info<< "    region " << regionI << " starts at " << cellI << endl;
+
+        // Make sure no parallel comms
+        bool oldParRun = UPstream::parRun();
+        UPstream::parRun() = false;
+
+        // Per region do a reordering.
+        fvMeshSubset subsetter(mesh);
+        subsetter.setLargeCellSubset(cellToRegion, regionI);
+
+        const fvMesh& subMesh = subsetter.subMesh();
+
+        labelList subReverseCellOrder = method.renumber
+        (
+            subMesh,
+            subMesh.cellCentres()
+        );
+
+        labelList subCellOrder
+        (
+            invert
+            (
+                subMesh.nCells(),
+                subReverseCellOrder
+            )
+        );
+
+        // Restore state
+        UPstream::parRun() = oldParRun;
+
+        const labelList& cellMap = subsetter.cellMap();
+
+        forAll(subCellOrder, i)
+        {
+            cellOrder[cellI++] = cellMap[subCellOrder[i]];
+        }
+    }
+    Info<< endl;
+
+    return cellOrder;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -505,23 +569,28 @@ int main(int argc, char *argv[])
     label blockSize = 0;
 
     // Construct renumberMethod
+    autoPtr<IOdictionary> renumberDictPtr;
     autoPtr<renumberMethod> renumberPtr;
 
     if (readDict)
     {
         Info<< "Renumber according to renumberMeshDict." << nl << endl;
 
-        IOdictionary renumberDict
+        renumberDictPtr.reset
         (
-            IOobject
+            new IOdictionary
             (
-                "renumberMeshDict",
-                runTime.system(),
-                mesh,
-                IOobject::MUST_READ_IF_MODIFIED,
-                IOobject::NO_WRITE
+                IOobject
+                (
+                    "renumberMeshDict",
+                    runTime.system(),
+                    mesh,
+                    IOobject::MUST_READ_IF_MODIFIED,
+                    IOobject::NO_WRITE
+                )
             )
         );
+        const IOdictionary renumberDict = renumberDictPtr();
 
         renumberPtr = renumberMethod::New(renumberDict);
 
@@ -562,7 +631,7 @@ int main(int argc, char *argv[])
                 << endl;
         }
 
-        writeMaps = readLabel(renumberDict.lookup("writeMaps"));
+        renumberDict.lookup("writeMaps") >> writeMaps;
         if (writeMaps)
         {
             Info<< "Writing renumber maps (new to old) to polyMesh." << nl
@@ -683,22 +752,15 @@ int main(int argc, char *argv[])
         // fields is done correctly!
 
         label nBlocks = mesh.nCells() / blockSize;
-        Info<< "nBlocks = " << nBlocks << endl;
+        Info<< "nBlocks   = " << nBlocks << endl;
 
-        // Read decomposePar dictionary
-        IOdictionary decomposeDict
-        (
-            IOobject
-            (
-                "decomposeParDict",
-                runTime.system(),
-                mesh,
-                IOobject::MUST_READ_IF_MODIFIED,
-                IOobject::NO_WRITE
-            )
-        );
+        // Read decompositionMethod dictionary
+        dictionary decomposeDict(renumberDictPtr().subDict("blockCoeffs"));
         decomposeDict.set("numberOfSubdomains", nBlocks);
 
+        bool oldParRun = UPstream::parRun();
+        UPstream::parRun() = false;
+
         autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
         (
             decomposeDict
@@ -713,6 +775,9 @@ int main(int argc, char *argv[])
             )
         );
 
+        // Restore state
+        UPstream::parRun() = oldParRun;
+
         // For debugging: write out region
         createScalarField
         (
@@ -726,23 +791,7 @@ int main(int argc, char *argv[])
             << nl << endl;
 
 
-        // Find point per region
-        pointField regionPoints(nBlocks, vector::zero);
-        forAll(cellToRegion, cellI)
-        {
-            regionPoints[cellToRegion[cellI]] = mesh.cellCentres()[cellI];
-        }
-
-        // Use block based renumbering.
-        // Detemines old to new cell ordering
-        labelList reverseCellOrder = renumberPtr().renumber
-        (
-            mesh,
-            cellToRegion,
-            regionPoints
-        );
-
-        cellOrder = invert(mesh.nCells(), reverseCellOrder);
+        cellOrder = regionRenumber(renumberPtr(), mesh, cellToRegion);
 
         // Determine new to old face order with new cell numbering
         faceOrder = getRegionFaceOrder
diff --git a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
index a852fad7db..988379720b 100644
--- a/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
+++ b/applications/utilities/mesh/manipulation/renumberMesh/renumberMeshDict
@@ -22,10 +22,14 @@ writeMaps true;
 // e.g. nonBlockingGaussSeidel.
 sortCoupledFaceCells false;
 
-// Optional entry: renumber on a block-by-block basis. This can be used on
-// large cases to keep the blocks fitting in cache with all the the cache
-// missed bunched at the end.
-//blockSize 0;
+// Optional entry: renumber on a block-by-block basis. It uses a
+// blockCoeffs dictionary to construct a decompositionMethod to do
+// a block subdivision) and then applies the renumberMethod to each
+// block in turn. This can be used in large cases to keep the blocks
+// fitting in cache with all the the cache misses bunched at the end.
+// This number is the approximate size of the blocks - this gets converted
+// to a number of blocks that is the input to the decomposition method.
+//blockSize 1000;
 
 // Optional entry: sort points into internal and boundary points
 //orderPoints false;
@@ -37,11 +41,11 @@ method          CuthillMcKee;
 //method          random;
 //method          spring;
 
-CuthillMcKeeCoeffs
-{
-    // Reverse CuthillMcKee (RCM) or plain
-    reverse true;
-}
+//CuthillMcKeeCoeffs
+//{
+//    // Reverse CuthillMcKee (RCM) or plain
+//    reverse true;
+//}
 
 
 manualCoeffs
@@ -65,4 +69,17 @@ springCoeffs
 }
 
 
+blockCoeffs
+{
+    method          scotch;
+    //method hierarchical;
+    //hierarchicalCoeffs
+    //{
+    //    n           (1 2 1);
+    //    delta       0.001;
+    //    order       xyz;
+    //}
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
index fd58570b34..02acd59236 100644
--- a/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
+++ b/applications/utilities/parallelProcessing/decomposePar/decomposePar.C
@@ -528,7 +528,7 @@ int main(int argc, char *argv[])
                         << "Cell number should be between 0 and "
                         << mesh.nCells()-1 << nl
                         << "On this mesh the particle should be in cell "
-                        <<  mesh.findCell(iter().position())
+                        << mesh.findCell(iter().position())
                         << exit(FatalError);
                 }
 
@@ -789,6 +789,14 @@ int main(int argc, char *argv[])
 
 
         // Point fields
+        if
+        (
+            pointScalarFields.size()
+         || pointVectorFields.size()
+         || pointSphericalTensorFields.size()
+         || pointSymmTensorFields.size()
+         || pointTensorFields.size()
+        )
         {
             labelIOList pointProcAddressing
             (
diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
index 5eaff912ef..cde7f5b12d 100644
--- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
+++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
@@ -795,9 +795,9 @@ int main(int argc, char *argv[])
     );
     fvMesh& mesh = meshPtr();
 
-    //Pout<< "Read mesh:" << endl;
-    //printMeshData(Pout, mesh);
-    //Pout<< endl;
+    // Print some statistics
+    Info<< "Before distribution:" << endl;
+    printMeshData(mesh);
 
 
 
@@ -1022,7 +1022,6 @@ int main(int argc, char *argv[])
     //map().distributeFaceData(faceCc);
 
 
-    // Print a bit
     // Print some statistics
     Info<< "After distribution:" << endl;
     printMeshData(mesh);
diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/Allwmake b/applications/utilities/postProcessing/graphics/PV3Readers/Allwmake
index 6851210bd7..1568f8afb8 100755
--- a/applications/utilities/postProcessing/graphics/PV3Readers/Allwmake
+++ b/applications/utilities/postProcessing/graphics/PV3Readers/Allwmake
@@ -1,6 +1,6 @@
 #!/bin/sh
 cd ${0%/*} || exit 1    # run from this directory
-set -x
+#set -x
 
 if [ -d "$ParaView_DIR" -a -r "$ParaView_DIR" ]
 then
@@ -15,6 +15,8 @@ then
     wmake libso vtkPV3Readers
     PV3blockMeshReader/Allwmake
     PV3FoamReader/Allwmake
+else
+    echo "ERROR: ParaView not found in $ParaView_DIR"
 fi
 
 # ----------------------------------------------------------------- end-of-file
diff --git a/applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_get_number_of_files_in_dataset.H b/applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_get_number_of_files_in_dataset.H
index e2246d27da..ed936ee082 100644
--- a/applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_get_number_of_files_in_dataset.H
+++ b/applications/utilities/postProcessing/graphics/ensightFoamReader/USERD_get_number_of_files_in_dataset.H
@@ -9,5 +9,5 @@ int USERD_get_number_of_files_in_dataset(void)
 
     // use 1 insted of 0 which gives an un-necessary warning.
     Num_dataset_files = 1;
-    return Num_dataset_files;;
+    return Num_dataset_files;
 }
diff --git a/applications/utilities/postProcessing/wall/wallShearStress/Make/options b/applications/utilities/postProcessing/wall/wallShearStress/Make/options
index 8862565863..b386fa4540 100644
--- a/applications/utilities/postProcessing/wall/wallShearStress/Make/options
+++ b/applications/utilities/postProcessing/wall/wallShearStress/Make/options
@@ -1,11 +1,14 @@
 EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/turbulenceModels \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    -lincompressibleRASModels \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels \
     -lfiniteVolume \
     -lgenericPatchFields
diff --git a/applications/utilities/postProcessing/wall/wallShearStress/createFields.H b/applications/utilities/postProcessing/wall/wallShearStress/createFields.H
deleted file mode 100644
index d069c4545d..0000000000
--- a/applications/utilities/postProcessing/wall/wallShearStress/createFields.H
+++ /dev/null
@@ -1,22 +0,0 @@
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhi.H"
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::RASModel> RASModel
-    (
-        incompressible::RASModel::New(U, phi, laminarTransport)
-    );
diff --git a/applications/utilities/postProcessing/wall/wallShearStress/wallShearStress.C b/applications/utilities/postProcessing/wall/wallShearStress/wallShearStress.C
index 7d440417df..6069837b95 100644
--- a/applications/utilities/postProcessing/wall/wallShearStress/wallShearStress.C
+++ b/applications/utilities/postProcessing/wall/wallShearStress/wallShearStress.C
@@ -25,23 +25,130 @@ Application
     wallShearStress
 
 Description
-    Calculates and writes the wall shear stress, for the specified times.
+    Calculates and reports wall shear stress for all patches, for the
+    specified times when using RAS turbulence models.
+
+    Default behaviour assumes operating in incompressible mode.
+    Use the -compressible option for compressible RAS cases.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+
 #include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "RASModel.H"
+#include "incompressible/RAS/RASModel/RASModel.H"
+
+#include "basicPsiThermo.H"
+#include "compressible/RAS/RASModel/RASModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+void calcIncompressible
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volVectorField& U,
+    volVectorField& wallShearStress
+)
+{
+    #include "createPhi.H"
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::RASModel> model
+    (
+        incompressible::RASModel::New(U, phi, laminarTransport)
+    );
+
+    const volSymmTensorField Reff(model->devReff());
+
+    forAll(wallShearStress.boundaryField(), patchI)
+    {
+        wallShearStress.boundaryField()[patchI] =
+        (
+           -mesh.Sf().boundaryField()[patchI]
+           /mesh.magSf().boundaryField()[patchI]
+        ) & Reff.boundaryField()[patchI];
+    }
+}
+
+
+void calcCompressible
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volVectorField& U,
+    volVectorField& wallShearStress
+)
+{
+    IOobject rhoHeader
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    );
+
+    if (!rhoHeader.headerOk())
+    {
+        Info<< "    no rho field" << endl;
+        return;
+    }
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho(rhoHeader, mesh);
+
+    #include "compressibleCreatePhi.H"
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    autoPtr<compressible::RASModel> model
+    (
+        compressible::RASModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    const volSymmTensorField Reff(model->devRhoReff());
+
+    forAll(wallShearStress.boundaryField(), patchI)
+    {
+        wallShearStress.boundaryField()[patchI] =
+        (
+           -mesh.Sf().boundaryField()[patchI]
+           /mesh.magSf().boundaryField()[patchI]
+        ) & Reff.boundaryField()[patchI];
+    }
+}
+
+
 int main(int argc, char *argv[])
 {
     timeSelector::addOptions();
+
+    #include "addRegionOption.H"
+
+    argList::addBoolOption
+    (
+        "compressible",
+        "calculate compressible wall shear stress"
+    );
+
     #include "setRootCase.H"
     #include "createTime.H"
     instantList timeDirs = timeSelector::select0(runTime, args);
-    #include "createMesh.H"
+    #include "createNamedMesh.H"
+
+    const bool compressible = args.optionFound("compressible");
 
     forAll(timeDirs, timeI)
     {
@@ -49,10 +156,6 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << endl;
         mesh.readUpdate();
 
-        #include "createFields.H"
-
-        volSymmTensorField Reff(RASModel->devReff());
-
         volVectorField wallShearStress
         (
             IOobject
@@ -67,19 +170,41 @@ int main(int argc, char *argv[])
             dimensionedVector
             (
                 "wallShearStress",
-                Reff.dimensions(),
+                sqr(dimLength)/sqr(dimTime),
                 vector::zero
             )
         );
 
-        forAll(wallShearStress.boundaryField(), patchi)
+        IOobject UHeader
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        );
+
+        if (UHeader.headerOk())
         {
-            wallShearStress.boundaryField()[patchi] =
-            (
-                -mesh.Sf().boundaryField()[patchi]
-                /mesh.magSf().boundaryField()[patchi]
-            ) & Reff.boundaryField()[patchi];
+            Info<< "Reading field U\n" << endl;
+            volVectorField U(UHeader, mesh);
+
+            if (compressible)
+            {
+                calcCompressible(mesh, runTime, U, wallShearStress);
+            }
+            else
+            {
+                calcIncompressible(mesh, runTime, U, wallShearStress);
+            }
         }
+        else
+        {
+            Info<< "    no U field" << endl;
+        }
+
+        Info<< "Writing wall shear stress to field " << wallShearStress.name()
+            << nl << endl;
 
         wallShearStress.write();
     }
diff --git a/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H b/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
index 2a203eac48..9015c473a2 100644
--- a/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
+++ b/applications/utilities/preProcessing/mapFields/MapConsistentVolFields.H
@@ -35,12 +35,13 @@ License
 namespace Foam
 {
 
-template<class Type>
+template<class Type, class CombineOp>
 void MapConsistentVolFields
 (
     const IOobjectList& objects,
     const meshToMesh& meshToMeshInterp,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const CombineOp& cop
 )
 {
     const fvMesh& meshSource = meshToMeshInterp.fromMesh();
@@ -84,7 +85,13 @@ void MapConsistentVolFields
             );
 
             // Interpolate field
-            meshToMeshInterp.interpolate(fieldTarget, fieldSource, mapOrder);
+            meshToMeshInterp.interpolate//<Type, eqOp<Type> >
+            (
+                fieldTarget,
+                fieldSource,
+                mapOrder,
+                cop
+            );
 
             // Write field
             fieldTarget.write();
@@ -97,7 +104,12 @@ void MapConsistentVolFields
             GeometricField<Type, fvPatchField, volMesh> fieldTarget
             (
                 fieldTargetIOobject,
-                meshToMeshInterp.interpolate(fieldSource, mapOrder)
+                meshToMeshInterp.interpolate//<Type, eqOp<Type> >
+                (
+                    fieldSource,
+                    mapOrder,
+                    cop
+                )
             );
 
             // Write field
diff --git a/applications/utilities/preProcessing/mapFields/MapMeshes.H b/applications/utilities/preProcessing/mapFields/MapMeshes.H
new file mode 100644
index 0000000000..e497eb0ccd
--- /dev/null
+++ b/applications/utilities/preProcessing/mapFields/MapMeshes.H
@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef MapMeshes_H
+#define MapMeshes_H
+
+#include "MapVolFields.H"
+#include "MapConsistentVolFields.H"
+#include "mapLagrangian.H"
+#include "UnMapped.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<template<class> class CombineOp>
+void MapConsistentMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh::order& mapOrder
+)
+{
+    // Create the interpolation scheme
+    meshToMesh meshToMeshInterp(meshSource, meshTarget);
+
+    Info<< nl
+        << "Consistently creating and mapping fields for time "
+        << meshSource.time().timeName() << nl << endl;
+
+    {
+        // Search for list of objects for this time
+        IOobjectList objects(meshSource, meshSource.time().timeName());
+
+        // Map volFields
+        // ~~~~~~~~~~~~~
+        MapConsistentVolFields<scalar>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<scalar>()
+        );
+        MapConsistentVolFields<vector>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<vector>()
+        );
+        MapConsistentVolFields<sphericalTensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<sphericalTensor>()
+        );
+        MapConsistentVolFields<symmTensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<symmTensor>()
+        );
+        MapConsistentVolFields<tensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<tensor>()
+        );
+    }
+
+    {
+        // Search for list of target objects for this time
+        IOobjectList objects(meshTarget, meshTarget.time().timeName());
+
+        // Mark surfaceFields as unmapped
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        UnMapped<surfaceScalarField>(objects);
+        UnMapped<surfaceVectorField>(objects);
+        UnMapped<surfaceSphericalTensorField>(objects);
+        UnMapped<surfaceSymmTensorField>(objects);
+        UnMapped<surfaceTensorField>(objects);
+
+        // Mark pointFields as unmapped
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        UnMapped<pointScalarField>(objects);
+        UnMapped<pointVectorField>(objects);
+        UnMapped<pointSphericalTensorField>(objects);
+        UnMapped<pointSymmTensorField>(objects);
+        UnMapped<pointTensorField>(objects);
+    }
+
+    mapLagrangian(meshToMeshInterp);
+}
+
+
+template<template<class> class CombineOp>
+void MapSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches,
+    const meshToMesh::order& mapOrder
+)
+{
+    // Create the interpolation scheme
+    meshToMesh meshToMeshInterp
+    (
+        meshSource,
+        meshTarget,
+        patchMap,
+        cuttingPatches
+    );
+
+    Info<< nl
+        << "Mapping fields for time " << meshSource.time().timeName()
+        << nl << endl;
+
+    {
+        // Search for list of source objects for this time
+        IOobjectList objects(meshSource, meshSource.time().timeName());
+
+        // Map volFields
+        // ~~~~~~~~~~~~~
+        MapVolFields<scalar>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<scalar>()
+        );
+        MapVolFields<vector>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<vector>()
+        );
+        MapVolFields<sphericalTensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<sphericalTensor>()
+        );
+        MapVolFields<symmTensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<symmTensor>()
+        );
+        MapVolFields<tensor>
+        (
+            objects,
+            meshToMeshInterp,
+            mapOrder,
+            CombineOp<tensor>()
+        );
+    }
+
+    {
+        // Search for list of target objects for this time
+        IOobjectList objects(meshTarget, meshTarget.time().timeName());
+
+        // Mark surfaceFields as unmapped
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        UnMapped<surfaceScalarField>(objects);
+        UnMapped<surfaceVectorField>(objects);
+        UnMapped<surfaceSphericalTensorField>(objects);
+        UnMapped<surfaceSymmTensorField>(objects);
+        UnMapped<surfaceTensorField>(objects);
+
+        // Mark pointFields as unmapped
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        UnMapped<pointScalarField>(objects);
+        UnMapped<pointVectorField>(objects);
+        UnMapped<pointSphericalTensorField>(objects);
+        UnMapped<pointSymmTensorField>(objects);
+        UnMapped<pointTensorField>(objects);
+    }
+
+    mapLagrangian(meshToMeshInterp);
+}
+
+
+template<template<class> class CombineOp>
+void MapConsistentSubMesh
+(
+    const fvMesh& meshSource,
+    const fvMesh& meshTarget,
+    const meshToMesh::order& mapOrder
+)
+{
+    HashTable<word> patchMap;
+    HashTable<label> cuttingPatchTable;
+
+    forAll(meshTarget.boundary(), patchi)
+    {
+        if (!isA<processorFvPatch>(meshTarget.boundary()[patchi]))
+        {
+            patchMap.insert
+            (
+                meshTarget.boundary()[patchi].name(),
+                meshTarget.boundary()[patchi].name()
+            );
+        }
+        else
+        {
+            cuttingPatchTable.insert
+            (
+                meshTarget.boundaryMesh()[patchi].name(),
+                -1
+            );
+        }
+    }
+
+    MapSubMesh<CombineOp>
+    (
+        meshSource,
+        meshTarget,
+        patchMap,
+        cuttingPatchTable.toc(),
+        mapOrder
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/preProcessing/mapFields/MapVolFields.H b/applications/utilities/preProcessing/mapFields/MapVolFields.H
index 64547a54de..65eb44d12d 100644
--- a/applications/utilities/preProcessing/mapFields/MapVolFields.H
+++ b/applications/utilities/preProcessing/mapFields/MapVolFields.H
@@ -35,12 +35,13 @@ License
 namespace Foam
 {
 
-template<class Type>
+template<class Type, class CombineOp>
 void MapVolFields
 (
     const IOobjectList& objects,
     const meshToMesh& meshToMeshInterp,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const CombineOp& cop
 )
 {
     const fvMesh& meshSource = meshToMeshInterp.fromMesh();
@@ -84,7 +85,13 @@ void MapVolFields
             );
 
             // Interpolate field
-            meshToMeshInterp.interpolate(fieldTarget, fieldSource, mapOrder);
+            meshToMeshInterp.interpolate//<Type, eqOp<Type> >
+            (
+                fieldTarget,
+                fieldSource,
+                mapOrder,
+                cop
+            );
 
             // Write field
             fieldTarget.write();
diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C
index 8d70911cab..9571be6fb0 100644
--- a/applications/utilities/preProcessing/mapFields/mapFields.C
+++ b/applications/utilities/preProcessing/mapFields/mapFields.C
@@ -34,11 +34,8 @@ Description
 
 #include "fvCFD.H"
 #include "meshToMesh.H"
-#include "MapVolFields.H"
-#include "MapConsistentVolFields.H"
-#include "UnMapped.H"
 #include "processorFvPatch.H"
-#include "mapLagrangian.H"
+#include "MapMeshes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -46,56 +43,28 @@ void mapConsistentMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const bool subtract
 )
 {
-    // Create the interpolation scheme
-    meshToMesh meshToMeshInterp(meshSource, meshTarget);
-
-    Info<< nl
-        << "Consistently creating and mapping fields for time "
-        << meshSource.time().timeName() << nl << endl;
-
+    if (subtract)
     {
-        // Search for list of objects for this time
-        IOobjectList objects(meshSource, meshSource.time().timeName());
-
-        // Map volFields
-        // ~~~~~~~~~~~~~
-        MapConsistentVolFields<scalar>(objects, meshToMeshInterp, mapOrder);
-        MapConsistentVolFields<vector>(objects, meshToMeshInterp, mapOrder);
-        MapConsistentVolFields<sphericalTensor>
+        MapConsistentMesh<minusEqOp>
         (
-            objects,
-            meshToMeshInterp,
+            meshSource,
+            meshTarget,
             mapOrder
         );
-        MapConsistentVolFields<symmTensor>(objects, meshToMeshInterp, mapOrder);
-        MapConsistentVolFields<tensor>(objects, meshToMeshInterp, mapOrder);
     }
-
+    else
     {
-        // Search for list of target objects for this time
-        IOobjectList objects(meshTarget, meshTarget.time().timeName());
-
-        // Mark surfaceFields as unmapped
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        UnMapped<surfaceScalarField>(objects);
-        UnMapped<surfaceVectorField>(objects);
-        UnMapped<surfaceSphericalTensorField>(objects);
-        UnMapped<surfaceSymmTensorField>(objects);
-        UnMapped<surfaceTensorField>(objects);
-
-        // Mark pointFields as unmapped
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        UnMapped<pointScalarField>(objects);
-        UnMapped<pointVectorField>(objects);
-        UnMapped<pointSphericalTensorField>(objects);
-        UnMapped<pointSymmTensorField>(objects);
-        UnMapped<pointTensorField>(objects);
+        MapConsistentMesh<eqOp>
+        (
+            meshSource,
+            meshTarget,
+            mapOrder
+        );
     }
-
-    mapLagrangian(meshToMeshInterp);
 }
 
 
@@ -105,57 +74,32 @@ void mapSubMesh
     const fvMesh& meshTarget,
     const HashTable<word>& patchMap,
     const wordList& cuttingPatches,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const bool subtract
 )
 {
-    // Create the interpolation scheme
-    meshToMesh meshToMeshInterp
-    (
-        meshSource,
-        meshTarget,
-        patchMap,
-        cuttingPatches
-    );
-
-    Info<< nl
-        << "Mapping fields for time " << meshSource.time().timeName()
-        << nl << endl;
-
+    if (subtract)
     {
-        // Search for list of source objects for this time
-        IOobjectList objects(meshSource, meshSource.time().timeName());
-
-        // Map volFields
-        // ~~~~~~~~~~~~~
-        MapVolFields<scalar>(objects, meshToMeshInterp, mapOrder);
-        MapVolFields<vector>(objects, meshToMeshInterp, mapOrder);
-        MapVolFields<sphericalTensor>(objects, meshToMeshInterp, mapOrder);
-        MapVolFields<symmTensor>(objects, meshToMeshInterp, mapOrder);
-        MapVolFields<tensor>(objects, meshToMeshInterp, mapOrder);
+        MapSubMesh<minusEqOp>
+        (
+            meshSource,
+            meshTarget,
+            patchMap,
+            cuttingPatches,
+            mapOrder
+        );
     }
-
+    else
     {
-        // Search for list of target objects for this time
-        IOobjectList objects(meshTarget, meshTarget.time().timeName());
-
-        // Mark surfaceFields as unmapped
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        UnMapped<surfaceScalarField>(objects);
-        UnMapped<surfaceVectorField>(objects);
-        UnMapped<surfaceSphericalTensorField>(objects);
-        UnMapped<surfaceSymmTensorField>(objects);
-        UnMapped<surfaceTensorField>(objects);
-
-        // Mark pointFields as unmapped
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        UnMapped<pointScalarField>(objects);
-        UnMapped<pointVectorField>(objects);
-        UnMapped<pointSphericalTensorField>(objects);
-        UnMapped<pointSymmTensorField>(objects);
-        UnMapped<pointTensorField>(objects);
+        MapSubMesh<eqOp>
+        (
+            meshSource,
+            meshTarget,
+            patchMap,
+            cuttingPatches,
+            mapOrder
+        );
     }
-
-    mapLagrangian(meshToMeshInterp);
 }
 
 
@@ -163,40 +107,28 @@ void mapConsistentSubMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
-    const meshToMesh::order& mapOrder
+    const meshToMesh::order& mapOrder,
+    const bool subtract
 )
 {
-    HashTable<word> patchMap;
-    HashTable<label> cuttingPatchTable;
-
-    forAll(meshTarget.boundary(), patchi)
+    if (subtract)
     {
-        if (!isA<processorFvPatch>(meshTarget.boundary()[patchi]))
-        {
-            patchMap.insert
-            (
-                meshTarget.boundary()[patchi].name(),
-                meshTarget.boundary()[patchi].name()
-            );
-        }
-        else
-        {
-            cuttingPatchTable.insert
-            (
-                meshTarget.boundaryMesh()[patchi].name(),
-                -1
-            );
-        }
+        MapConsistentSubMesh<minusEqOp>
+        (
+            meshSource,
+            meshTarget,
+            mapOrder
+        );
+    }
+    else
+    {
+        MapConsistentSubMesh<eqOp>
+        (
+            meshSource,
+            meshTarget,
+            mapOrder
+        );
     }
-
-    mapSubMesh
-    (
-        meshSource,
-        meshTarget,
-        patchMap,
-        cuttingPatchTable.toc(),
-        mapOrder
-    );
 }
 
 
@@ -288,6 +220,11 @@ int main(int argc, char *argv[])
         "word",
         "specify the mapping method"
     );
+    argList::addBoolOption
+    (
+        "subtract",
+        "subtract mapped source from target"
+    );
 
     argList args(argc, argv);
 
@@ -350,6 +287,13 @@ int main(int argc, char *argv[])
         Info<< "Mapping method: " << mapMethod << endl;
     }
 
+    const bool subtract = args.optionFound("subtract");
+    if (subtract)
+    {
+        Info<< "Subtracting mapped source field from target" << endl;
+    }
+
+
     #include "createTimes.H"
 
     HashTable<word> patchMap;
@@ -431,7 +375,13 @@ int main(int argc, char *argv[])
 
             if (consistent)
             {
-                mapConsistentSubMesh(meshSource, meshTarget, mapOrder);
+                mapConsistentSubMesh
+                (
+                    meshSource,
+                    meshTarget,
+                    mapOrder,
+                    subtract
+                );
             }
             else
             {
@@ -441,7 +391,8 @@ int main(int argc, char *argv[])
                     meshTarget,
                     patchMap,
                     cuttingPatches,
-                    mapOrder
+                    mapOrder,
+                    subtract
                 );
             }
         }
@@ -503,7 +454,13 @@ int main(int argc, char *argv[])
 
             if (consistent)
             {
-                mapConsistentSubMesh(meshSource, meshTarget, mapOrder);
+                mapConsistentSubMesh
+                (
+                    meshSource,
+                    meshTarget,
+                    mapOrder,
+                    subtract
+                );
             }
             else
             {
@@ -513,7 +470,8 @@ int main(int argc, char *argv[])
                     meshTarget,
                     patchMap,
                     addProcessorPatches(meshTarget, cuttingPatches),
-                    mapOrder
+                    mapOrder,
+                    subtract
                 );
             }
         }
@@ -629,7 +587,8 @@ int main(int argc, char *argv[])
                             (
                                 meshSource,
                                 meshTarget,
-                                mapOrder
+                                mapOrder,
+                                subtract
                             );
                         }
                         else
@@ -640,7 +599,8 @@ int main(int argc, char *argv[])
                                 meshTarget,
                                 patchMap,
                                 addProcessorPatches(meshTarget, cuttingPatches),
-                                mapOrder
+                                mapOrder,
+                                subtract
                             );
                         }
                     }
@@ -679,7 +639,7 @@ int main(int argc, char *argv[])
 
         if (consistent)
         {
-            mapConsistentMesh(meshSource, meshTarget, mapOrder);
+            mapConsistentMesh(meshSource, meshTarget, mapOrder, subtract);
         }
         else
         {
@@ -689,7 +649,8 @@ int main(int argc, char *argv[])
                 meshTarget,
                 patchMap,
                 cuttingPatches,
-                mapOrder
+                mapOrder,
+                subtract
             );
         }
     }
diff --git a/src/Allwmake b/src/Allwmake
index a0622b30f7..97023aa700 100755
--- a/src/Allwmake
+++ b/src/Allwmake
@@ -46,7 +46,7 @@ wmake $makeType genericPatchFields
 # Build the proper scotchDecomp, metisDecomp etc.
 parallel/Allwmake $*
 
-wmake $makeType renumberMethods
+renumber/Allwmake $*
 
 wmake $makeType conversion
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
index afd31e642a..76f82fa78c 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H
@@ -144,11 +144,13 @@ public:
             inline tetPointRef oldTet(const polyMesh& mesh) const;
 
             //- Return the geometry corresponding to the tri on the
-            //  mesh face for this tet from the supplied mesh
+            //  mesh face for this tet from the supplied mesh. Normal of
+            //  the tri points out of the cell.
             inline triPointRef faceTri(const polyMesh& mesh) const;
 
             //- Return the point indices corresponding to the tri on the mesh
-            //  face for this tet from the supplied mesh
+            //  face for this tet from the supplied mesh. Normal of
+            //  the tri points out of the cell.
             inline triFace faceTriIs(const polyMesh& mesh) const;
 
             //- Return the geometry corresponding to the tri on the
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index 3867915943..496442012d 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
@@ -277,9 +277,9 @@ Foam::scalar Foam::triangle<Point, PointRef>::barycentric
 
     bary.setSize(3);
 
-    bary[0] = (d11*d20 - d01*d21)/denom;
-    bary[1] = (d00*d21 - d01*d20)/denom;
-    bary[2] = 1.0 - bary[0] - bary[1];
+    bary[1] = (d11*d20 - d01*d21)/denom;
+    bary[2] = (d00*d21 - d01*d20)/denom;
+    bary[0] = 1.0 - bary[1] - bary[2];
 
     return denom;
 }
diff --git a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
index 98a742e16b..7ff31db0ec 100644
--- a/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
+++ b/src/finiteVolume/cfdTools/general/solutionControl/pimpleControl/pimpleControl.C
@@ -196,9 +196,10 @@ bool Foam::pimpleControl::loop()
     bool completed = false;
     if (criteriaSatisfied())
     {
-        Info<< algorithmName_ << ": converged in " << corr_ << " iterations"
+        Info<< algorithmName_ << ": converged in " << corr_ - 1 << " iterations"
             << endl;
         completed = true;
+        corr_ = 0;
     }
     else
     {
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
index 3ef207b702..13828ff7cf 100644
--- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
+++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/cellPointWeight/cellPointWeight.C
@@ -83,8 +83,8 @@ void Foam::cellPointWeight::findTetrahedron
             )
             {
                 faceVertices_[0] = f[tetIs.faceBasePt()];
-                faceVertices_[1] = f[tetIs.facePtA()];;
-                faceVertices_[2] = f[tetIs.facePtB()];;
+                faceVertices_[1] = f[tetIs.facePtA()];
+                faceVertices_[2] = f[tetIs.facePtB()];
 
                 return;
             }
@@ -191,8 +191,8 @@ void Foam::cellPointWeight::findTriangle
                 weights_[3] = triWeights[2];
 
                 faceVertices_[0] = f[tetIs.faceBasePt()];
-                faceVertices_[1] = f[tetIs.facePtA()];;
-                faceVertices_[2] = f[tetIs.facePtB()];;
+                faceVertices_[1] = f[tetIs.facePtA()];
+                faceVertices_[2] = f[tetIs.facePtB()];
 
                 return;
             }
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
index 296029b7d5..48eff18670 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -286,6 +286,10 @@ bool Foam::meshRefinement::isCollapsedFace
     const label faceI
 ) const
 {
+    // Severe nonorthogonality threshold
+    const scalar severeNonorthogonalityThreshold =
+        ::cos(degToRad(maxNonOrtho));
+
     vector s = mesh_.faces()[faceI].normal(points);
     scalar magS = mag(s);
 
@@ -301,11 +305,11 @@ bool Foam::meshRefinement::isCollapsedFace
     if (mesh_.isInternalFace(faceI))
     {
         label nei = mesh_.faceNeighbour()[faceI];
-        vector d = ownCc - mesh_.cellCentres()[nei];
+        vector d = mesh_.cellCentres()[nei] - ownCc;
 
         scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
 
-        if (dDotS < maxNonOrtho)
+        if (dDotS < severeNonorthogonalityThreshold)
         {
             return true;
         }
@@ -320,11 +324,11 @@ bool Foam::meshRefinement::isCollapsedFace
 
         if (mesh_.boundaryMesh()[patchI].coupled())
         {
-            vector d = ownCc - neiCc[faceI-mesh_.nInternalFaces()];
+            vector d = neiCc[faceI-mesh_.nInternalFaces()] - ownCc;
 
             scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
 
-            if (dDotS < maxNonOrtho)
+            if (dDotS < severeNonorthogonalityThreshold)
             {
                 return true;
             }
@@ -337,7 +341,7 @@ bool Foam::meshRefinement::isCollapsedFace
         {
             // Collapsing normal boundary face does not cause problems with
             // non-orthogonality
-            return true;
+            return false;
         }
     }
 }
@@ -663,7 +667,8 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                 )
                 {
                     nPrevented++;
-                    //Pout<< "Preventing collapse of 8 anchor point cell "
+                    //Pout<< "Preventing baffling/removal of 8 anchor point"
+                    //    << " cell "
                     //    << cellI << " at " << mesh_.cellCentres()[cellI]
                     //    << " since new volume "
                     //    << mesh_.cells()[cellI].mag(newPoints, mesh_.faces())
@@ -747,7 +752,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                     )
                     {
                         nPrevented++;
-                        //Pout<< "Preventing collapse of 7 anchor cell "
+                        //Pout<< "Preventing baffling of 7 anchor cell "
                         //    << cellI
                         //    << " at " << mesh_.cellCentres()[cellI]
                         //    << " since new volume "
@@ -853,7 +858,8 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                 )
                 {
                     nPrevented++;
-                    //Pout<< "Preventing collapse of face " << faceI
+                    //Pout<< "Preventing baffling (to avoid collapse) of face "
+                    //    << faceI
                     //    << " with all boundary edges "
                     //    << " at " << mesh_.faceCentres()[faceI]
                     //    << endl;
@@ -909,7 +915,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
                         )
                         {
                             nPrevented++;
-                            //Pout<< "Preventing collapse of coupled face "
+                            //Pout<< "Preventing baffling of coupled face "
                             //    << faceI
                             //    << " with all boundary edges "
                             //    << " at " << mesh_.faceCentres()[faceI]
@@ -943,7 +949,7 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     {
         Info<< "markFacesOnProblemCells : prevented "
             << returnReduce(nPrevented, sumOp<label>())
-            << " internal faces fom getting converted into baffles."
+            << " internal faces from getting converted into baffles."
             << endl;
     }
 
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index 363000bad6..e0e96412d3 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -236,8 +236,8 @@ void Foam::decompositionMethod::calcCellCells
                 label globalNei = globalNeighbour[bFaceI];
                 if
                 (
-                    globalAgglom.isLocal(globalNei)
-                 && globalAgglom.toLocal(globalNei) != own
+                   !globalAgglom.isLocal(globalNei)
+                 || globalAgglom.toLocal(globalNei) != own
                 )
                 {
                     nFacesPerCell[own]++;
@@ -285,10 +285,11 @@ void Foam::decompositionMethod::calcCellCells
                 label own = agglom[faceOwner[faceI]];
 
                 label globalNei = globalNeighbour[bFaceI];
+
                 if
                 (
-                    globalAgglom.isLocal(globalNei)
-                 && globalAgglom.toLocal(globalNei) != own
+                   !globalAgglom.isLocal(globalNei)
+                 || globalAgglom.toLocal(globalNei) != own
                 )
                 {
                     m[offsets[own] + nFacesPerCell[own]++] = globalNei;
@@ -304,15 +305,15 @@ void Foam::decompositionMethod::calcCellCells
     // Check for duplicates connections between cells
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Done as postprocessing step since we now have cellCells.
-    label startIndex = 0;
     label newIndex = 0;
     labelHashSet nbrCells;
     forAll(cellCells, cellI)
     {
         nbrCells.clear();
-        nbrCells.insert(cellI);
+        nbrCells.insert(globalAgglom.toGlobal(cellI));
 
-        label& endIndex = cellCells.offsets()[cellI+1];
+        label startIndex = cellCells.offsets()[cellI];
+        label endIndex = cellCells.offsets()[cellI+1];
 
         for (label i = startIndex; i < endIndex; i++)
         {
@@ -321,9 +322,7 @@ void Foam::decompositionMethod::calcCellCells
                 cellCells.m()[newIndex++] = cellCells.m()[i];
             }
         }
-
-        startIndex = endIndex;
-        endIndex = newIndex;
+        cellCells.offsets()[cellI+1] = newIndex;
     }
 
     cellCells.m().setSize(newIndex);
diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
index 54be6ad834..e7c4995a36 100644
--- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
+++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C
@@ -371,29 +371,6 @@ Foam::label Foam::ptscotchDecomp::decompose
         Pout<< "ptscotchDecomp : entering with xadj:" << xadj.size() << endl;
     }
 
-
-if (debug)
-{
-    Pout<< "nProcessors_:" << nProcessors_ << endl;
-
-    globalIndex globalCells(xadj.size()-1);
-
-    Pout<< "Xadj:" << endl;
-    for (label cellI = 0; cellI < xadj.size()-1; cellI++)
-    {
-        Pout<< "cell:" << cellI
-            << "  global:" << globalCells.toGlobal(cellI)
-            << " connected to:" << endl;
-        label start = xadj[cellI];
-        label end = xadj[cellI+1];
-        for (label i = start; i < end; i++)
-        {
-            Pout<< "    cell:" << adjncy[i] << endl;
-        }
-    }
-    Pout<< endl;
-}
-
     // Dump graph
     if (decompositionDict_.found("ptscotchCoeffs"))
     {
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.C b/src/postProcessing/functionObjects/field/streamLine/streamLine.C
index 30d8052513..ab815a8df8 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLine.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.C
@@ -34,6 +34,7 @@ License
 #include "globalIndex.H"
 #include "mapDistribute.H"
 #include "interpolationCellPoint.H"
+#include "PatchTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -41,6 +42,57 @@ defineTypeNameAndDebug(Foam::streamLine, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+Foam::autoPtr<Foam::indirectPrimitivePatch>
+Foam::streamLine::wallPatch() const
+{
+    const fvMesh& mesh = dynamic_cast<const fvMesh&>(obr_);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    label nFaces = 0;
+
+    forAll(patches, patchI)
+    {
+        //if (!polyPatch::constraintType(patches[patchI].type()))
+        if (isA<wallPolyPatch>(patches[patchI]))
+        {
+            nFaces += patches[patchI].size();
+        }
+    }
+
+    labelList addressing(nFaces);
+
+    nFaces = 0;
+
+    forAll(patches, patchI)
+    {
+        //if (!polyPatch::constraintType(patches[patchI].type()))
+        if (isA<wallPolyPatch>(patches[patchI]))
+        {
+            const polyPatch& pp = patches[patchI];
+
+            forAll(pp, i)
+            {
+                addressing[nFaces++] = pp.start()+i;
+            }
+        }
+    }
+
+    return autoPtr<indirectPrimitivePatch>
+    (
+        new indirectPrimitivePatch
+        (
+            IndirectList<face>
+            (
+                mesh.faces(),
+                addressing
+            ),
+            mesh.points()
+        )
+    );
+}
+
+
 void Foam::streamLine::track()
 {
     const Time& runTime = obr_.time();
@@ -72,7 +124,7 @@ void Foam::streamLine::track()
 
     label nSeeds = returnReduce(particles.size(), sumOp<label>());
 
-    Info<< "Seeded " << nSeeds << " particles." << endl;
+    Info<< type() << " : seeded " << nSeeds << " particles." << endl;
 
     // Read or lookup fields
     PtrList<volScalarField> vsFlds;
@@ -163,7 +215,6 @@ void Foam::streamLine::track()
                 vsInterp.set
                 (
                     nScalar++,
-                    //new interpolationCellPoint<scalar>(f)
                     interpolation<scalar>::New
                     (
                         interpolationScheme_,
@@ -186,7 +237,6 @@ void Foam::streamLine::track()
                 vvInterp.set
                 (
                     nVector++,
-                    //new interpolationCellPoint<vector>(f)
                     interpolation<vector>::New
                     (
                         interpolationScheme_,
@@ -241,6 +291,120 @@ void Foam::streamLine::track()
         allVectors_[i].setCapacity(nSeeds);
     }
 
+
+    //- For surface following: per face the normal to constrain the velocity
+    //  with.
+    pointField patchNormals;
+
+
+    if (followSurface_)
+    {
+        // Determine geometric data on the non-constraint patches
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        // - point normals
+        // - neighbouring cells
+        autoPtr<indirectPrimitivePatch> boundaryPatch(wallPatch());
+
+        // Calculate parallel consistent normal (on boundaryPatch points only)
+        pointField boundaryNormals
+        (
+            PatchTools::pointNormals
+            (
+                mesh,
+                boundaryPatch(),
+                boundaryPatch().addressing()
+            )
+        );
+        // Make sure all points know about point normal. There might be points
+        // which are not on boundaryPatch() but coupled to points that are.
+        pointField pointNormals(mesh.nPoints(), vector::zero);
+        UIndirectList<point>
+        (
+            pointNormals,
+            boundaryPatch().meshPoints()
+        ) = boundaryNormals;
+        boundaryNormals.clear();
+
+        mesh.globalData().syncPointData
+        (
+            pointNormals,
+            maxMagSqrEqOp<point>(),
+            mapDistribute::transform()
+        );
+
+        // Calculate per-face a single constraint normal:
+        // - faces on patch: face-normal
+        // - faces on cell with point on patch: point-normal
+        // - rest: vector::zero
+        //
+        // Note: the disadvantage is that we can only handle faces with one
+        //       point on the patch - otherwise we might pick up the wrong
+        //       velocity.
+
+
+
+        patchNormals.setSize(mesh.nFaces(), vector::zero);
+        label nPointNormals = 0;
+
+        // 1. faces on patch
+        UIndirectList<point>(patchNormals, boundaryPatch().addressing()) =
+             boundaryPatch().faceNormals();
+
+        // 2. faces with a point on the patch
+        forAll(mesh.faces(), faceI)
+        {
+            if (patchNormals[faceI] == vector::zero)
+            {
+                const face& f = mesh.faces()[faceI];
+                forAll(f, fp)
+                {
+                    label pointI = f[fp];
+                    if (pointNormals[pointI] != vector::zero)
+                    {
+                        patchNormals[faceI] = pointNormals[pointI];
+                        nPointNormals++;
+                    }
+                }
+            }
+        }
+
+        // 3. cells on faces with a point on the patch
+        forAll(mesh.points(), pointI)
+        {
+            if (pointNormals[pointI] != vector::zero)
+            {
+                const labelList& pCells = mesh.pointCells(pointI);
+                forAll(pCells, i)
+                {
+                    const cell& cFaces = mesh.cells()[pCells[i]];
+                    forAll(cFaces, j)
+                    {
+                        label faceI = cFaces[j];
+
+                        if (patchNormals[faceI] == vector::zero)
+                        {
+                            patchNormals[faceI] = pointNormals[pointI];
+                            nPointNormals++;
+                        }
+                    }
+                }
+            }
+        }
+
+        Info<< type() << " : calculated constrained-tracking normals" << nl
+            << "Number of faces in mesh                     : "
+            << returnReduce(mesh.nFaces(), sumOp<label>()) << nl
+            << "    of which on walls                       : "
+            <<  returnReduce
+                (
+                    boundaryPatch().addressing().size(),
+                    sumOp<label>()
+                ) << nl
+            << "    of which on cell with point on boundary : "
+            << returnReduce(nPointNormals, sumOp<label>()) << nl
+            << endl;
+    }
+
     // additional particle info
     streamLineParticle::trackingData td
     (
@@ -248,8 +412,16 @@ void Foam::streamLine::track()
         vsInterp,
         vvInterp,
         UIndex,         // index of U in vvInterp
-        trackForward_,  // track in +u direction
-        nSubCycle_,     // step through cells in steps?
+        trackForward_,  // track in +u direction?
+        nSubCycle_,     // automatic track control:step through cells in steps?
+        trackLength_,   // fixed track length
+
+        //- For surface-constrained streamlines
+        followSurface_, // stay on surface?
+        patchNormals,   // per face normal to constrain tracking to
+        //pointNormals,   // per point          ,,
+        trackHeight_,   // track height
+
         allTracks_,
         allScalars_,
         allVectors_
@@ -279,7 +451,8 @@ Foam::streamLine::streamLine
     name_(name),
     obr_(obr),
     loadFromFiles_(loadFromFiles),
-    active_(true)
+    active_(true),
+    nSubCycle_(0)
 {
     // Only active if a fvMesh is available
     if (isA<fvMesh>(obr_))
@@ -334,6 +507,16 @@ void Foam::streamLine::read(const dictionary& dict)
                 dict.lookup("U") >> UName_;
             }
         }
+
+        if (findIndex(fields_, UName_) == -1)
+        {
+            FatalIOErrorIn("streamLine::read(const dictionary&)", dict)
+                << "Velocity field for tracking " << UName_
+                << " should be present in the list of fields " << fields_
+                << exit(FatalIOError);
+        }
+
+
         dict.lookup("trackForward") >> trackForward_;
         dict.lookup("lifeTime") >> lifeTime_;
         if (lifeTime_ < 1)
@@ -343,10 +526,48 @@ void Foam::streamLine::read(const dictionary& dict)
                 << exit(FatalError);
         }
 
-        dict.lookup("nSubCycle") >> nSubCycle_;
-        if (nSubCycle_ < 1)
+
+        bool subCycling = dict.found("nSubCycle");
+        bool fixedLength = dict.found("trackLength");
+
+        if (subCycling && fixedLength)
         {
-            nSubCycle_ = 1;
+            FatalIOErrorIn("streamLine::read(const dictionary&)", dict)
+                << "Cannot both specify automatic time stepping (through '"
+                << "nSubCycle' specification) and fixed track length (through '"
+                << "trackLength')"
+                << exit(FatalIOError);
+        }
+
+
+        nSubCycle_ = 1;
+        if (dict.readIfPresent("nSubCycle", nSubCycle_))
+        {
+            trackLength_ = VGREAT;
+            if (nSubCycle_ < 1)
+            {
+                nSubCycle_ = 1;
+            }
+            Info<< type() << " : automatic track length specified through"
+                << " number of sub cycles : " << nSubCycle_ << nl << endl;
+        }
+        else
+        {
+            dict.lookup("trackLength") >> trackLength_;
+
+            Info<< type() << " : fixed track length specified : "
+                << trackLength_ << nl << endl;
+        }
+
+
+        followSurface_ = false;
+        trackHeight_ = VGREAT;
+        if (dict.readIfPresent("followSurface", followSurface_))
+        {
+            dict.lookup("trackHeight") >> trackHeight_;
+
+            Info<< type() << " : surface-constrained streamline at "
+                << trackHeight_ << " distance above the surface" << nl << endl;
         }
 
         interpolationScheme_ = dict.lookupOrDefault
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLine.H b/src/postProcessing/functionObjects/field/streamLine/streamLine.H
index 9fc981034c..77984ed2c8 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLine.H
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLine.H
@@ -43,6 +43,7 @@ SourceFiles
 #include "vectorList.H"
 #include "polyMesh.H"
 #include "writer.H"
+#include "indirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -95,9 +96,17 @@ class streamLine
         //- Maximum lifetime (= number of cells) of particle
         label lifeTime_;
 
+        //- Surface-constrained?
+        bool followSurface_;
+
+        scalar trackHeight_;
+
         //- Number of subcycling steps
         label nSubCycle_;
 
+        //- Track length
+        scalar trackLength_;
+
         //- Optional specified name of particles
         word cloudName_;
 
@@ -141,6 +150,8 @@ class streamLine
             List<DynamicList<vectorList> > allVectors_;
 
 
+        //- Construct patch out of all wall patch faces
+        autoPtr<indirectPrimitivePatch> wallPatch() const;
 
         //- Do all seeding and tracking
         void track();
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
index c7808b9900..279e17ef80 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.C
@@ -43,7 +43,7 @@ Foam::scalar Foam::streamLineParticle::calcSubCycleDeltaT
     const vector& U
 ) const
 {
-    streamLineParticle testParticle(*this);
+    particle testParticle(*this);
 
     bool oldKeepParticle = td.keepParticle;
     bool oldSwitchProcessor = td.switchProcessor;
@@ -59,7 +59,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
 (
     const trackingData& td,
     const point& position,
-    const label cellI
+    const label cellI,
+    const label faceI
 )
 {
     if (cellI == -1)
@@ -76,7 +77,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
             td.vsInterp_[scalarI].interpolate
             (
                 position,
-                cellI
+                cellI,
+                faceI
             )
         );
     }
@@ -89,7 +91,8 @@ Foam::vector Foam::streamLineParticle::interpolateFields
             td.vvInterp_[vectorI].interpolate
             (
                 position,
-                cellI
+                cellI,
+                faceI
             )
         );
     }
@@ -201,13 +204,24 @@ bool Foam::streamLineParticle::move
 
             // Store current position and sampled velocity.
             sampledPositions_.append(position());
-            vector U = interpolateFields(td, position(), cell());
+            vector U = interpolateFields(td, position(), cell(), face());
 
             if (!td.trackForward_)
             {
                 U = -U;
             }
 
+            if (td.followSurface_)
+            {
+                //- Simple: remove wall-normal component. Should keep particle
+                //  at same height. Does not work well on warped faces.
+                const vector& n = td.patchNormals_[tetFaceI_];
+                if (n != vector::zero)
+                {
+                    U -= (U&n)*n;
+                }
+            }
+
             scalar magU = mag(U);
 
             if (magU < SMALL)
@@ -219,7 +233,13 @@ bool Foam::streamLineParticle::move
 
             U /= magU;
 
-            if (subIter == 0 && td.nSubCycle_ > 1)
+            if (td.trackLength_ < GREAT)
+            {
+                dt = td.trackLength_;
+                //Pout<< "    subiteration " << subIter
+                //    << " : fixed length: updated dt:" << dt << endl;
+            }
+            else if (subIter == 0 && td.nSubCycle_ > 1)
             {
                 // Adapt dt to cross cell in a few steps
                 dt = calcSubCycleDeltaT(td, dt, U);
@@ -264,7 +284,8 @@ bool Foam::streamLineParticle::move
             if (debug)
             {
                 Pout<< "streamLineParticle : Removing stagnant particle:"
-                    << p << " sampled positions:" << sampledPositions_.size()
+                    << p.position()
+                    << " sampled positions:" << sampledPositions_.size()
                     << endl;
             }
             td.keepParticle = false;
@@ -273,12 +294,13 @@ bool Foam::streamLineParticle::move
         {
             // Normal exit. Store last position and fields
             sampledPositions_.append(position());
-            interpolateFields(td, position(), cell());
+            interpolateFields(td, position(), cell(), face());
 
             if (debug)
             {
                 Pout<< "streamLineParticle : Removing particle:"
-                    << p << " sampled positions:" << sampledPositions_.size()
+                    << p.position()
+                    << " sampled positions:" << sampledPositions_.size()
                     << endl;
             }
         }
@@ -374,8 +396,42 @@ void Foam::streamLineParticle::hitWallPatch
     const tetIndices&
 )
 {
-    // Remove particle
-    td.keepParticle = false;
+    if (td.followSurface_)
+    {
+        // Push slightly in. Since we're pushing to tet centre we should still
+        // be in the same tet.
+        tetPointRef tet = currentTet();
+
+        const point oldPosition(position());
+
+        position() += trackingCorrectionTol*(tet.centre() - position());
+
+        if (debug)
+        {
+            label patchID = patch(face());
+
+            Pout<< "hitWallPatch : on wall "
+                << mesh().boundaryMesh()[patchID].name()
+                << " moved from " << oldPosition
+                << " to " << position()
+                << " due to centre " << tet.centre() << endl;
+        }
+
+        // Check that still in tet (should always be the case
+        if (debug && !tet.inside(position()))
+        {
+            FatalErrorIn("streamLineParticle::hitWallPatch()")
+                << "Pushing particle " << *this
+                << " away from wall " << wpp.name()
+                << " moved it outside of tet."
+                << exit(FatalError);
+        }
+    }
+    else
+    {
+        // Remove particle
+        td.keepParticle = false;
+    }
 }
 
 
diff --git a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
index b1ee9acf88..c5feff4285 100644
--- a/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
+++ b/src/postProcessing/functionObjects/field/streamLine/streamLineParticle.H
@@ -73,6 +73,13 @@ public:
         const label UIndex_;
         const bool trackForward_;
         const label nSubCycle_;
+        const scalar trackLength_;
+
+        //- For surface-constrained tracking
+        const bool followSurface_;
+        const pointField& patchNormals_;
+        const scalar trackHeight_;
+
 
         DynamicList<vectorList>& allPositions_;
         List<DynamicList<scalarList> >& allScalars_;
@@ -89,6 +96,12 @@ public:
                 const label UIndex,
                 const bool trackForward,
                 const label nSubCycle,
+                const scalar trackLength,
+                const bool followSurface,
+
+                const pointField& patchNormals,
+                const scalar trackHeight,
+
                 DynamicList<List<point> >& allPositions,
                 List<DynamicList<scalarList> >& allScalars,
                 List<DynamicList<vectorList> >& allVectors
@@ -100,6 +113,12 @@ public:
                 UIndex_(UIndex),
                 trackForward_(trackForward),
                 nSubCycle_(nSubCycle),
+                trackLength_(trackLength),
+                followSurface_(followSurface),
+
+                patchNormals_(patchNormals),
+                trackHeight_(trackHeight),
+
                 allPositions_(allPositions),
                 allScalars_(allScalars),
                 allVectors_(allVectors)
@@ -135,12 +154,20 @@ private:
             const vector& U
         ) const;
 
+        void constrainVelocity
+        (
+            trackingData& td,
+            const scalar dt,
+            vector& U
+        );
+
         //- Interpolate all quantities; return interpolated velocity.
         vector interpolateFields
         (
             const trackingData&,
             const point&,
-            const label cellI
+            const label cellI,
+            const label faceI
         );
 
 
diff --git a/src/renumber/Allwmake b/src/renumber/Allwmake
new file mode 100755
index 0000000000..36c057dafe
--- /dev/null
+++ b/src/renumber/Allwmake
@@ -0,0 +1,30 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+makeType=${1:-libso}
+
+## get ZOLTAN_ARCH_PATH
+#if settings=`$WM_PROJECT_DIR/bin/foamEtcFile config/zoltan.sh`
+#then
+#    . $settings
+#    echo "using ZOLTAN_ARCH_PATH=$ZOLTAN_ARCH_PATH"
+#else
+#    echo
+#    echo "Error: no config/zoltan.sh settings"
+#    echo
+#fi
+
+
+set -x
+
+wmake $makeType renumberMethods
+
+#if [ -n "$ZOLTAN_ARCH_PATH" ]
+#then
+#    wmake $makeType zoltanRenumber
+#else
+#    echo
+#    echo "Skipping zoltanRenumber"
+#    echo
+#fi
+
+# ----------------------------------------------------------------- end-of-file
diff --git a/src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
similarity index 90%
rename from src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
rename to src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
index e1ec010335..2732484b36 100644
--- a/src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
+++ b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
@@ -63,7 +63,7 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 (
     const polyMesh& mesh,
     const pointField& points
-)
+) const
 {
     CompactListList<label> cellCells;
     decompositionMethod::calcCellCells
@@ -79,10 +79,7 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 
     if (reverse_)
     {
-        for (label i = 0; i < orderedToOld.size()/2; i++)
-        {
-            Swap(orderedToOld[i], orderedToOld[orderedToOld.size()-i-1]);
-        }
+        reverse(orderedToOld);
     }
 
     return invert(orderedToOld.size(), orderedToOld);
@@ -93,16 +90,13 @@ Foam::labelList Foam::CuthillMcKeeRenumber::renumber
 (
     const labelListList& cellCells,
     const pointField& points
-)
+) const
 {
     labelList orderedToOld = bandCompression(cellCells);
 
     if (reverse_)
     {
-        for (label i = 0; i < orderedToOld.size()/2; i++)
-        {
-            Swap(orderedToOld[i], orderedToOld[orderedToOld.size()-i-1]);
-        }
+        reverse(orderedToOld);
     }
 
     return invert(orderedToOld.size(), orderedToOld);
diff --git a/src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
similarity index 97%
rename from src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
rename to src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
index 03148c8739..a4ec011b40 100644
--- a/src/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
+++ b/src/renumber/renumberMethods/CuthillMcKeeRenumber/CuthillMcKeeRenumber.H
@@ -82,7 +82,7 @@ public:
 
         //- Return for every coordinate the wanted processor number.
         //  We need a polyMesh (to be able to load the file)
-        virtual labelList renumber(const pointField&)
+        virtual labelList renumber(const pointField&) const
         {
             notImplemented("CuthillMcKeeRenumber::renumber(const pointField&)");
             return labelList(0);
@@ -94,7 +94,7 @@ public:
         (
             const polyMesh& mesh,
             const pointField& cc
-        );
+        ) const;
 
         //- Return for every cell the new cell label.
         //  The connectivity is equal to mesh.cellCells() except
@@ -103,7 +103,7 @@ public:
         (
             const labelListList& cellCells,
             const pointField& cc
-        );
+        ) const;
 };
 
 
diff --git a/src/renumberMethods/Make/files b/src/renumber/renumberMethods/Make/files
similarity index 82%
rename from src/renumberMethods/Make/files
rename to src/renumber/renumberMethods/Make/files
index fb366f73c4..c19f1e6b6d 100644
--- a/src/renumberMethods/Make/files
+++ b/src/renumber/renumberMethods/Make/files
@@ -3,5 +3,6 @@ manualRenumber/manualRenumber.C
 CuthillMcKeeRenumber/CuthillMcKeeRenumber.C
 randomRenumber/randomRenumber.C
 springRenumber/springRenumber.C
+boundaryFirstRenumber/boundaryFirstRenumber.C
 
 LIB = $(FOAM_LIBBIN)/librenumberMethods
diff --git a/src/renumberMethods/Make/options b/src/renumber/renumberMethods/Make/options
similarity index 100%
rename from src/renumberMethods/Make/options
rename to src/renumber/renumberMethods/Make/options
diff --git a/src/renumber/renumberMethods/boundaryFirstRenumber/boundaryFirstRenumber.C b/src/renumber/renumberMethods/boundaryFirstRenumber/boundaryFirstRenumber.C
new file mode 100644
index 0000000000..293e5ba19b
--- /dev/null
+++ b/src/renumber/renumberMethods/boundaryFirstRenumber/boundaryFirstRenumber.C
@@ -0,0 +1,201 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "boundaryFirstRenumber.H"
+#include "addToRunTimeSelectionTable.H"
+#include "bandCompression.H"
+#include "decompositionMethod.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(boundaryFirstRenumber, 0);
+
+    addToRunTimeSelectionTable
+    (
+        renumberMethod,
+        boundaryFirstRenumber,
+        dictionary
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::boundaryFirstRenumber::boundaryFirstRenumber
+(
+    const dictionary& renumberDict
+)
+:
+    renumberMethod(renumberDict),
+    reverse_
+    (
+        renumberDict.found(typeName + "Coeffs")
+      ? Switch(renumberDict.subDict(typeName + "Coeffs").lookup("reverse"))
+      : Switch(true)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::labelList Foam::boundaryFirstRenumber::renumber
+(
+    const polyMesh& mesh,
+    const pointField& points
+) const
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+
+    // Distance to boundary. Minimum of neighbours.
+    labelList distance(mesh.nCells(), -1);
+
+    // New order
+    labelList newOrder(mesh.nCells());
+    label cellInOrder = 0;
+
+
+    // Starting faces for walk. These are the zero distance faces.
+    DynamicList<label> frontFaces(mesh.nFaces()-mesh.nInternalFaces());
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        //- Note: cannot check for empty since these are introduced by
+        //  fvMeshSubset. Also most loops don't care about patch type.
+        //if (!isA<emptyPolyPatch>(pp))
+        {
+            forAll(pp, i)
+            {
+                frontFaces.append(pp.start()+i);
+            }
+        }
+    }
+
+    // Loop over all frontFaces
+    label currentDistance = 0;
+    while (frontFaces.size() > 0)
+    {
+        DynamicList<label> frontCells(frontFaces.size());
+
+        // Set all of frontFaces' neighbours to current distance
+
+        forAll(frontFaces, i)
+        {
+            label faceI = frontFaces[i];
+
+            if (mesh.isInternalFace(faceI))
+            {
+                label ownCellI = own[faceI];
+                if (distance[ownCellI] == -1)
+                {
+                    distance[ownCellI] = currentDistance;
+                    frontCells.append(ownCellI);
+                }
+                label neiCellI = nei[faceI];
+                if (distance[neiCellI] == -1)
+                {
+                    distance[neiCellI] = currentDistance;
+                    frontCells.append(neiCellI);
+                }
+            }
+            else
+            {
+                label ownCellI = own[faceI];
+                if (distance[ownCellI] == -1)
+                {
+                    distance[ownCellI] = currentDistance;
+                    frontCells.append(ownCellI);
+                }
+            }
+        }
+
+        //Pout<< "For distance:" << currentDistance
+        //    << " from " << frontFaces.size() << " faces to "
+        //    << frontCells.size() << " cells." << endl;
+
+
+        // TBD. Determine order within current shell (frontCells). For now
+        // just add them.
+        forAll(frontCells, i)
+        {
+            newOrder[cellInOrder] = frontCells[i];
+            cellInOrder++;
+        }
+
+        // From cells to faces
+        frontFaces.clear();
+        forAll(frontCells, i)
+        {
+            label cellI = frontCells[i];
+            const cell& cFaces = mesh.cells()[cellI];
+
+            forAll(cFaces, i)
+            {
+                label faceI = cFaces[i];
+                if (mesh.isInternalFace(faceI))
+                {
+                    label nbrCellI =
+                    (
+                        mesh.faceOwner()[faceI] == cellI
+                      ? mesh.faceNeighbour()[faceI]
+                      : mesh.faceOwner()[faceI]
+                    );
+                    if (distance[nbrCellI] == -1)
+                    {
+                        frontFaces.append(faceI);
+                    }
+                }
+            }
+        }
+
+        currentDistance++;
+    }
+
+    // Return furthest away cell first
+    if (reverse_)
+    {
+        reverse(newOrder);
+    }
+
+    //forAll(newOrder, i)
+    //{
+    //    label cellI = newOrder[i];
+    //
+    //    Pout<< "cell:" << cellI << endl;
+    //    Pout<< " at distance:" << distance[cellI]
+    //        << endl;
+    //}
+
+    return invert(newOrder.size(), newOrder);
+}
+
+
+// ************************************************************************* //
diff --git a/src/renumber/renumberMethods/boundaryFirstRenumber/boundaryFirstRenumber.H b/src/renumber/renumberMethods/boundaryFirstRenumber/boundaryFirstRenumber.H
new file mode 100644
index 0000000000..690af21455
--- /dev/null
+++ b/src/renumber/renumberMethods/boundaryFirstRenumber/boundaryFirstRenumber.H
@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::boundaryFirstRenumber
+
+Description
+    Cuthill-McKee renumbering
+
+SourceFiles
+    boundaryFirstRenumber.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef boundaryFirstRenumber_H
+#define boundaryFirstRenumber_H
+
+#include "renumberMethod.H"
+#include "Switch.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class boundaryFirstRenumber Declaration
+\*---------------------------------------------------------------------------*/
+
+class boundaryFirstRenumber
+:
+    public renumberMethod
+{
+    // Private data
+
+        const Switch reverse_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        void operator=(const boundaryFirstRenumber&);
+        boundaryFirstRenumber(const boundaryFirstRenumber&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("boundaryFirst");
+
+
+    // Constructors
+
+        //- Construct given the renumber dictionary
+        boundaryFirstRenumber(const dictionary& renumberDict);
+
+
+    //- Destructor
+    virtual ~boundaryFirstRenumber()
+    {}
+
+
+    // Member Functions
+
+        //- Return for every coordinate the wanted processor number.
+        //  We need a polyMesh (to be able to load the file)
+        virtual labelList renumber(const pointField&) const
+        {
+            notImplemented
+            (
+                "boundaryFirstRenumber::renumber(const pointField&)"
+            );
+            return labelList(0);
+        }
+
+        //- Return for every coordinate the wanted processor number. Use the
+        //  mesh connectivity (if needed)
+        virtual labelList renumber
+        (
+            const polyMesh& mesh,
+            const pointField& cc
+        ) const;
+
+        //- Return for every cell the new cell label.
+        //  The connectivity is equal to mesh.cellCells() except
+        //  - the connections are across coupled patches
+        virtual labelList renumber
+        (
+            const labelListList& cellCells,
+            const pointField& cc
+        ) const
+        {
+            notImplemented
+            (
+                "boundaryFirstRenumber::renumber"
+                "(const labelListList&, const pointField&)"
+            );
+            return labelList(0);
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/renumberMethods/manualRenumber/manualRenumber.C b/src/renumber/renumberMethods/manualRenumber/manualRenumber.C
similarity index 99%
rename from src/renumberMethods/manualRenumber/manualRenumber.C
rename to src/renumber/renumberMethods/manualRenumber/manualRenumber.C
index 40989264b5..92b5c8e00c 100644
--- a/src/renumberMethods/manualRenumber/manualRenumber.C
+++ b/src/renumber/renumberMethods/manualRenumber/manualRenumber.C
@@ -61,7 +61,7 @@ Foam::labelList Foam::manualRenumber::renumber
 (
     const polyMesh& mesh,
     const pointField& points
-)
+) const
 {
     labelIOList oldToNew
     (
diff --git a/src/renumberMethods/manualRenumber/manualRenumber.H b/src/renumber/renumberMethods/manualRenumber/manualRenumber.H
similarity index 97%
rename from src/renumberMethods/manualRenumber/manualRenumber.H
rename to src/renumber/renumberMethods/manualRenumber/manualRenumber.H
index 635e93c614..a4d4b2f46f 100644
--- a/src/renumberMethods/manualRenumber/manualRenumber.H
+++ b/src/renumber/renumberMethods/manualRenumber/manualRenumber.H
@@ -81,7 +81,7 @@ public:
 
         //- Return for every coordinate the wanted processor number.
         //  We need a polyMesh (to be able to load the file)
-        virtual labelList renumber(const pointField&)
+        virtual labelList renumber(const pointField&) const
         {
             notImplemented("manualRenumber::renumber(const pointField&)");
             return labelList(0);
@@ -93,7 +93,7 @@ public:
         (
             const polyMesh& mesh,
             const pointField& cc
-        );
+        ) const;
 
         //- Return for every cell the new cell label.
         //  The connectivity is equal to mesh.cellCells() except
@@ -102,7 +102,7 @@ public:
         (
             const labelListList& cellCells,
             const pointField& cc
-        )
+        ) const
         {
             notImplemented
             (
diff --git a/src/renumberMethods/randomRenumber/randomRenumber.C b/src/renumber/renumberMethods/randomRenumber/randomRenumber.C
similarity index 99%
rename from src/renumberMethods/randomRenumber/randomRenumber.C
rename to src/renumber/renumberMethods/randomRenumber/randomRenumber.C
index 27ddc29c3f..2dae28e438 100644
--- a/src/renumberMethods/randomRenumber/randomRenumber.C
+++ b/src/renumber/renumberMethods/randomRenumber/randomRenumber.C
@@ -55,7 +55,7 @@ Foam::randomRenumber::randomRenumber(const dictionary& renumberDict)
 Foam::labelList Foam::randomRenumber::renumber
 (
     const pointField& points
-)
+) const
 {
     Random rndGen(0);
 
@@ -77,7 +77,7 @@ Foam::labelList Foam::randomRenumber::renumber
 (
     const polyMesh& mesh,
     const pointField& points
-)
+) const
 {
     return renumber(points);
 }
@@ -87,7 +87,7 @@ Foam::labelList Foam::randomRenumber::renumber
 (
     const labelListList& cellCells,
     const pointField& points
-)
+) const
 {
     return renumber(points);
 }
diff --git a/src/renumberMethods/randomRenumber/randomRenumber.H b/src/renumber/renumberMethods/randomRenumber/randomRenumber.H
similarity index 97%
rename from src/renumberMethods/randomRenumber/randomRenumber.H
rename to src/renumber/renumberMethods/randomRenumber/randomRenumber.H
index be704d1097..0665d5159b 100644
--- a/src/renumberMethods/randomRenumber/randomRenumber.H
+++ b/src/renumber/renumberMethods/randomRenumber/randomRenumber.H
@@ -76,7 +76,7 @@ public:
 
         //- Return for every coordinate the wanted processor number.
         //  We need a polyMesh (to be able to load the file)
-        virtual labelList renumber(const pointField&);
+        virtual labelList renumber(const pointField&) const;
 
         //- Return for every coordinate the wanted processor number. Use the
         //  mesh connectivity (if needed)
@@ -84,7 +84,7 @@ public:
         (
             const polyMesh& mesh,
             const pointField& cc
-        );
+        ) const;
 
         //- Return for every cell the new cell label.
         //  The connectivity is equal to mesh.cellCells() except
@@ -93,7 +93,7 @@ public:
         (
             const labelListList& cellCells,
             const pointField& cc
-        );
+        ) const;
 };
 
 
diff --git a/src/renumberMethods/renumberMethod/renumberMethod.C b/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
similarity index 99%
rename from src/renumberMethods/renumberMethod/renumberMethod.C
rename to src/renumber/renumberMethods/renumberMethod/renumberMethod.C
index 53138673b8..76ba1911e2 100644
--- a/src/renumberMethods/renumberMethod/renumberMethod.C
+++ b/src/renumber/renumberMethods/renumberMethod/renumberMethod.C
@@ -72,7 +72,7 @@ Foam::labelList Foam::renumberMethod::renumber
 (
     const polyMesh& mesh,
     const pointField& points
-)
+) const
 {
     CompactListList<label> cellCells;
     decompositionMethod::calcCellCells
@@ -94,7 +94,7 @@ Foam::labelList Foam::renumberMethod::renumber
     const polyMesh& mesh,
     const labelList& fineToCoarse,
     const pointField& coarsePoints
-)
+) const
 {
     CompactListList<label> coarseCellCells;
     decompositionMethod::calcCellCells
diff --git a/src/renumberMethods/renumberMethod/renumberMethod.H b/src/renumber/renumberMethods/renumberMethod/renumberMethod.H
similarity index 97%
rename from src/renumberMethods/renumberMethod/renumberMethod.H
rename to src/renumber/renumberMethods/renumberMethod/renumberMethod.H
index b4b6e89f22..ead92b1e41 100644
--- a/src/renumberMethods/renumberMethod/renumberMethod.H
+++ b/src/renumber/renumberMethods/renumberMethod/renumberMethod.H
@@ -111,7 +111,7 @@ public:
 
         //- Return for every cell the new cell label.
         //  This is only defined for geometric renumberMethods.
-        virtual labelList renumber(const pointField&)
+        virtual labelList renumber(const pointField&) const
         {
             notImplemented
             (
@@ -122,7 +122,7 @@ public:
 
         //- Return for every cell the new cell label. Use the
         //  mesh connectivity (if needed)
-        virtual labelList renumber(const polyMesh&, const pointField&);
+        virtual labelList renumber(const polyMesh&, const pointField&) const;
 
         //- Return for every cell the new cell label. Gets
         //  passed agglomeration map (from fine to coarse cells) and coarse
@@ -136,7 +136,7 @@ public:
             const polyMesh& mesh,
             const labelList& cellToRegion,
             const pointField& regionPoints
-        );
+        ) const;
 
         //- Return for every cell the new cell label.
         //  The connectivity is equal to mesh.cellCells() except
@@ -145,7 +145,7 @@ public:
         (
             const labelListList& cellCells,
             const pointField& cc
-        ) = 0;
+        ) const = 0;
 
 };
 
diff --git a/src/renumberMethods/springRenumber/springRenumber.C b/src/renumber/renumberMethods/springRenumber/springRenumber.C
similarity index 99%
rename from src/renumberMethods/springRenumber/springRenumber.C
rename to src/renumber/renumberMethods/springRenumber/springRenumber.C
index 5c6fa49350..64a15145ce 100644
--- a/src/renumberMethods/springRenumber/springRenumber.C
+++ b/src/renumber/renumberMethods/springRenumber/springRenumber.C
@@ -60,7 +60,7 @@ Foam::labelList Foam::springRenumber::renumber
 (
     const polyMesh& mesh,
     const pointField& points
-)
+) const
 {
     CompactListList<label> cellCells;
     decompositionMethod::calcCellCells
@@ -80,7 +80,7 @@ Foam::labelList Foam::springRenumber::renumber
 (
     const labelListList& cellCells,
     const pointField& points
-)
+) const
 {
     // Look at cell index as a 1D position parameter.
     // Move cells to the average 'position' of their neighbour.
diff --git a/src/renumberMethods/springRenumber/springRenumber.H b/src/renumber/renumberMethods/springRenumber/springRenumber.H
similarity index 97%
rename from src/renumberMethods/springRenumber/springRenumber.H
rename to src/renumber/renumberMethods/springRenumber/springRenumber.H
index 64342af1ec..30f7a3d838 100644
--- a/src/renumberMethods/springRenumber/springRenumber.H
+++ b/src/renumber/renumberMethods/springRenumber/springRenumber.H
@@ -96,7 +96,7 @@ public:
     // Member Functions
 
         //- Return for every coordinate the wanted processor number.
-        virtual labelList renumber(const pointField&)
+        virtual labelList renumber(const pointField&) const
         {
             notImplemented("springRenumber::renumber(const pointField&)");
             return labelList(0);
@@ -108,7 +108,7 @@ public:
         (
             const polyMesh& mesh,
             const pointField& cc
-        );
+        ) const;
 
         //- Return for every cell the new cell label.
         //  The connectivity is equal to mesh.cellCells() except
@@ -117,7 +117,7 @@ public:
         (
             const labelListList& cellCells,
             const pointField& cc
-        );
+        ) const;
 };
 
 
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
index 785829cc60..83d1ee902a 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
@@ -219,84 +219,93 @@ public:
         // Interpolation
 
             //- Map field
-            template<class Type>
+            template<class Type, class CombineOp>
             void mapField
             (
                 Field<Type>&,
                 const Field<Type>&,
-                const labelList& adr
+                const labelList& adr,
+                const CombineOp& cop
             ) const;
 
             //- Interpolate field using inverse-distance weights
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolateField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
                 const labelList& adr,
-                const scalarListList& weights
+                const scalarListList& weights,
+                const CombineOp& cop
             ) const;
 
             //- Interpolate field using cell-point interpolation
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolateField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
                 const labelList& adr,
-                const vectorField& centres
+                const vectorField& centres,
+                const CombineOp& cop
             ) const;
 
 
             //- Interpolate internal volume field
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolateInternalField
             (
                 Field<Type>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolateInternalField
             (
                 Field<Type>&,
                 const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
 
             //- Interpolate volume field
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolate
             (
                 GeometricField<Type, fvPatchField, volMesh>&,
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
-            template<class Type>
+            template<class Type, class CombineOp>
             void interpolate
             (
                 GeometricField<Type, fvPatchField, volMesh>&,
                 const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
 
             //- Interpolate volume field
-            template<class Type>
+            template<class Type, class CombineOp>
             tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
             (
                 const GeometricField<Type, fvPatchField, volMesh>&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 
-            template<class Type>
+            template<class Type, class CombineOp>
             tmp<GeometricField<Type, fvPatchField, volMesh> > interpolate
             (
                 const tmp<GeometricField<Type, fvPatchField, volMesh> >&,
-                order=INTERPOLATE
+                order=INTERPOLATE,
+                const CombineOp& cop = eqOp<Type>()
             ) const;
 };
 
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
index 991ae31826..7ed543ed3c 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C
@@ -31,12 +31,13 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::mapField
 (
     Field<Type>& toF,
     const Field<Type>& fromVf,
-    const labelList& adr
+    const labelList& adr,
+    const CombineOp& cop
 ) const
 {
     // Direct mapping of nearest-cell values
@@ -45,7 +46,7 @@ void Foam::meshToMesh::mapField
     {
         if (adr[celli] != -1)
         {
-            toF[celli] = fromVf[adr[celli]];
+            cop(toF[celli], fromVf[adr[celli]]);
         }
     }
 
@@ -53,13 +54,14 @@ void Foam::meshToMesh::mapField
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
     const labelList& adr,
-    const scalarListList& weights
+    const scalarListList& weights,
+    const CombineOp& cop
 ) const
 {
     // Inverse distance weighted interpolation
@@ -74,24 +76,27 @@ void Foam::meshToMesh::interpolateField
             const labelList& neighbours = cc[adr[celli]];
             const scalarList& w = weights[celli];
 
-            toF[celli] = fromVf[adr[celli]]*w[0];
+            Type f = fromVf[adr[celli]]*w[0];
 
             for (label ni = 1; ni < w.size(); ni++)
             {
-                toF[celli] += fromVf[neighbours[ni - 1]]*w[ni];
+                f += fromVf[neighbours[ni - 1]]*w[ni];
             }
+
+            cop(toF[celli], f);
         }
     }
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
     const labelList& adr,
-    const vectorField& centres
+    const vectorField& centres,
+    const CombineOp& cop
 ) const
 {
     // Cell-Point interpolation
@@ -101,31 +106,36 @@ void Foam::meshToMesh::interpolateField
     {
         if (adr[celli] != -1)
         {
-            toF[celli] = interpolator.interpolate
+            cop
             (
-                centres[celli],
-                adr[celli]
+                toF[celli],
+                interpolator.interpolate
+                (
+                    centres[celli],
+                    adr[celli]
+                )
             );
         }
     }
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateInternalField
 (
     Field<Type>& toF,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
     if (fromVf.mesh() != fromMesh_)
     {
         FatalErrorIn
         (
-            "meshToMesh::interpolateInternalField(Field<Type>& toF, "
-            "const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
-            "meshToMesh::order ord) const"
+            "meshToMesh::interpolateInternalField(Field<Type>&, "
+            "const GeometricField<Type, fvPatchField, volMesh>&, "
+            "meshToMesh::order, const CombineOp&) const"
         )   << "the argument field does not correspond to the right mesh. "
             << "Field size: " << fromVf.size()
             << " mesh size: " << fromMesh_.nCells()
@@ -136,9 +146,9 @@ void Foam::meshToMesh::interpolateInternalField
     {
         FatalErrorIn
         (
-            "meshToMesh::interpolateInternalField(Field<Type>& toF, "
-            "const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
-            "meshToMesh::order ord) const"
+            "meshToMesh::interpolateInternalField(Field<Type>&, "
+            "const GeometricField<Type, fvPatchField, volMesh>&, "
+            "meshToMesh::order, const CombineOp&) const"
         )   << "the argument field does not correspond to the right mesh. "
             << "Field size: " << toF.size()
             << " mesh size: " << toMesh_.nCells()
@@ -148,7 +158,7 @@ void Foam::meshToMesh::interpolateInternalField
     switch(ord)
     {
         case MAP:
-            mapField(toF, fromVf, cellAddressing_);
+            mapField(toF, fromVf, cellAddressing_, cop);
         break;
 
         case INTERPOLATE:
@@ -157,7 +167,8 @@ void Foam::meshToMesh::interpolateInternalField
                 toF,
                 fromVf,
                 cellAddressing_,
-                inverseDistanceWeights()
+                inverseDistanceWeights(),
+                cop
             );
         break;
 
@@ -167,44 +178,47 @@ void Foam::meshToMesh::interpolateInternalField
                 toF,
                 fromVf,
                 cellAddressing_,
-                toMesh_.cellCentres()
+                toMesh_.cellCentres(),
+                cop
             );
         break;
 
         default:
             FatalErrorIn
             (
-                "meshToMesh::interpolateInternalField(Field<Type>& toF, "
-                "const GeometricField<Type, fvPatchField, volMesh>& fromVf, "
-                "meshToMesh::order ord) const"
+                "meshToMesh::interpolateInternalField(Field<Type>&, "
+                "const GeometricField<Type, fvPatchField, volMesh>&, "
+                "meshToMesh::order, const CombineOp&) const"
             )   << "unknown interpolation scheme " << ord
                 << exit(FatalError);
     }
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateInternalField
 (
     Field<Type>& toF,
     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
-    interpolateInternalField(toF, tfromVf(), ord);
+    interpolateInternalField(toF, tfromVf(), ord, cop);
     tfromVf.clear();
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolate
 (
     GeometricField<Type, fvPatchField, volMesh>& toVf,
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
-    interpolateInternalField(toVf, fromVf, ord);
+    interpolateInternalField(toVf, fromVf, ord, cop);
 
     forAll(toMesh_.boundaryMesh(), patchi)
     {
@@ -219,7 +233,8 @@ void Foam::meshToMesh::interpolate
                     (
                         toVf.boundaryField()[patchi],
                         fromVf,
-                        boundaryAddressing_[patchi]
+                        boundaryAddressing_[patchi],
+                        cop
                     );
                 break;
 
@@ -229,7 +244,8 @@ void Foam::meshToMesh::interpolate
                         toVf.boundaryField()[patchi],
                         fromVf,
                         boundaryAddressing_[patchi],
-                        toPatch.Cf()
+                        toPatch.Cf(),
+                        cop
                     );
                 break;
 
@@ -239,7 +255,8 @@ void Foam::meshToMesh::interpolate
                         toVf.boundaryField()[patchi],
                         fromVf,
                         boundaryAddressing_[patchi],
-                        toPatch.Cf()
+                        toPatch.Cf(),
+                        cop
                     );
                 break;
 
@@ -247,9 +264,9 @@ void Foam::meshToMesh::interpolate
                     FatalErrorIn
                     (
                         "meshToMesh::interpolate("
-                        "GeometricField<Type, fvPatchField, volMesh>& toVf, "
-                        "const GeometricField<Type, fvPatchField, volMesh>& "
-                        "fromVf, meshToMesh::order ord) const"
+                        "GeometricField<Type, fvPatchField, volMesh>&, "
+                        "const GeometricField<Type, fvPatchField, volMesh>&, "
+                        "meshToMesh::order, const CombineOp&) const"
                     )   << "unknown interpolation scheme " << ord
                         << exit(FatalError);
             }
@@ -286,37 +303,40 @@ void Foam::meshToMesh::interpolate
                 [
                     fromMeshPatches_.find(patchMap_.find(toPatch.name())())()
                 ],
-                boundaryAddressing_[patchi]
+                boundaryAddressing_[patchi],
+                cop
             );
         }
     }
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolate
 (
     GeometricField<Type, fvPatchField, volMesh>& toVf,
     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
-    interpolate(toVf, tfromVf(), ord);
+    interpolate(toVf, tfromVf(), ord, cop);
     tfromVf.clear();
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 Foam::tmp< Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 Foam::meshToMesh::interpolate
 (
     const GeometricField<Type, fvPatchField, volMesh>& fromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
     // Create and map the internal-field values
     Field<Type> internalField(toMesh_.nCells());
-    interpolateInternalField(internalField, fromVf, ord);
+    interpolateInternalField(internalField, fromVf, ord, cop);
 
     // check whether both meshes have got the same number
     // of boundary patches
@@ -325,8 +345,8 @@ Foam::meshToMesh::interpolate
         FatalErrorIn
         (
             "meshToMesh::interpolate"
-            "(const GeometricField<Type, fvPatchField, volMesh>& fromVf,"
-            "meshToMesh::order ord) const"
+            "(const GeometricField<Type, fvPatchField, volMesh>&,"
+            "meshToMesh::order, const CombineOp&) const"
         )   << "Incompatible meshes: different number of boundaries, "
                "only internal field may be interpolated"
             << exit(FatalError);
@@ -381,16 +401,17 @@ Foam::meshToMesh::interpolate
 }
 
 
-template<class Type>
+template<class Type, class CombineOp>
 Foam::tmp< Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
 Foam::meshToMesh::interpolate
 (
     const tmp<GeometricField<Type, fvPatchField, volMesh> >& tfromVf,
-    meshToMesh::order ord
+    meshToMesh::order ord,
+    const CombineOp& cop
 ) const
 {
     tmp<GeometricField<Type, fvPatchField, volMesh> > tint =
-        interpolate(tfromVf(), ord);
+        interpolate(tfromVf(), ord, cop);
     tfromVf.clear();
 
     return tint;
diff --git a/tutorials/incompressible/simpleFoam/motorBike/system/streamLines b/tutorials/incompressible/simpleFoam/motorBike/system/streamLines
index 8b82eccbb1..50b9eefdd1 100644
--- a/tutorials/incompressible/simpleFoam/motorBike/system/streamLines
+++ b/tutorials/incompressible/simpleFoam/motorBike/system/streamLines
@@ -28,8 +28,15 @@ streamLines
     // Steps particles can travel before being removed
     lifeTime        10000;
 
-    // Number of steps per cell (estimate). Set to 1 to disable subcycling.
-    nSubCycle 5;
+    //- Specify either absolute length of steps (trackLength) or a number
+    //  of subcycling steps per cell (nSubCycle)
+
+        // Size of single track segment [m]
+        //trackLength 1e-3;
+
+        // Number of steps per cell (estimate). Set to 1 to disable subcycling.
+        nSubCycle 5;
+
 
     // Cloud name to use
     cloudName       particleTracks;
diff --git a/tutorials/incompressible/simpleFoam/turbineSiting/system/snappyHexMeshDict b/tutorials/incompressible/simpleFoam/turbineSiting/system/snappyHexMeshDict
index 1e6f477fc4..8054858f2a 100644
--- a/tutorials/incompressible/simpleFoam/turbineSiting/system/snappyHexMeshDict
+++ b/tutorials/incompressible/simpleFoam/turbineSiting/system/snappyHexMeshDict
@@ -173,7 +173,7 @@ castellatedMeshControls
         }
     }
 
-    resolveFeatureAngle 2;
+    resolveFeatureAngle 30;
 
 
     // Region-wise refinement