diff --git a/applications/solvers/combustion/XiFoam/XiEngineFoam/startSummary.H b/applications/solvers/combustion/XiFoam/XiEngineFoam/startSummary.H
index 2d0c6fd707..f8755e15cc 100644
--- a/applications/solvers/combustion/XiFoam/XiEngineFoam/startSummary.H
+++ b/applications/solvers/combustion/XiFoam/XiEngineFoam/startSummary.H
@@ -8,7 +8,7 @@ if (Pstream::master())
     (
         new OFstream
         (
-            runTime.rootPath()/runTime.globalCaseName()
+            runTime.globalPath()
            /("logSummary." + runTime.timeName() + ".dat")
         )
     );
diff --git a/applications/solvers/combustion/coldEngineFoam/startSummary.H b/applications/solvers/combustion/coldEngineFoam/startSummary.H
index b4b72be696..93c8d9493f 100644
--- a/applications/solvers/combustion/coldEngineFoam/startSummary.H
+++ b/applications/solvers/combustion/coldEngineFoam/startSummary.H
@@ -8,7 +8,7 @@ if (Pstream::master())
     (
         new OFstream
         (
-            runTime.rootPath()/runTime.globalCaseName()
+            runTime.globalPath()
            /("logSummary." + runTime.timeName() + ".dat")
         )
     );
diff --git a/applications/test/globalMeshData/Test-globalMeshData.C b/applications/test/globalMeshData/Test-globalMeshData.C
index 56a285ae02..8d3f2c24e6 100644
--- a/applications/test/globalMeshData/Test-globalMeshData.C
+++ b/applications/test/globalMeshData/Test-globalMeshData.C
@@ -180,7 +180,7 @@ int main(int argc, char *argv[])
         const labelListList& transformedSlaves =
             globalData.globalPointTransformedBoundaryFaces();
 
-        label nBnd = mesh.nFaces()-mesh.nInternalFaces();
+        const label nBnd = mesh.nBoundaryFaces();
 
         pointField fc(globalPointBoundaryFacesMap.constructSize());
         SubList<point>(fc, nBnd) =
diff --git a/applications/test/syncTools/Test-syncTools.C b/applications/test/syncTools/Test-syncTools.C
index 5dfef2fcf2..5fca906370 100644
--- a/applications/test/syncTools/Test-syncTools.C
+++ b/applications/test/syncTools/Test-syncTools.C
@@ -195,7 +195,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         SubList<face>
         (
             mesh.faces(),
-            mesh.nFaces()-mesh.nInternalFaces(),
+            mesh.nBoundaryFaces(),
             mesh.nInternalFaces()
         ),
         mesh.points()
diff --git a/applications/test/tensor/Test-tensor.C b/applications/test/tensor/Test-tensor.C
index d597985da5..cd72200975 100644
--- a/applications/test/tensor/Test-tensor.C
+++ b/applications/test/tensor/Test-tensor.C
@@ -1,4 +1,5 @@
 #include "tensor.H"
+#include "triad.H"
 #include "symmTensor.H"
 #include "transform.H"
 #include "stringList.H"
@@ -13,7 +14,40 @@ int main()
 
     tensor t3 = t1 + t2;
 
-    Info<< t3 << endl;
+    Info<< "tensor " << t1 << " + " << t2 << nl
+        << " = " << t3 << nl << nl;
+
+    {
+        Info<< "rows:" << nl;
+        for (direction i=0; i < 3; ++i)
+        {
+            Info<< "  (" << i << ") = " << t3.row(i) << nl;
+        }
+    }
+    {
+        Info<< "cols:" << nl;
+        for (direction i=0; i < 3; ++i)
+        {
+            Info<< "  (" << i << ") = " << t3.col(i) << nl;
+        }
+        Info<< "col<0> = " << t3.col<0>() << nl;
+        Info<< "col<1> = " << t3.col<1>() << nl;
+        Info<< "col<2> = " << t3.col<2>() << nl;
+        // Compilation error:  Info << "col<3> = " << t3.col<3>() << nl;
+
+        t3.col<1>({0, 2, 4});
+        Info<< "replaced col<1> = " << t3.col<1>() << nl;
+        Info<< "tensor " << t3 << nl;
+
+        t3.row<2>(Zero);
+        Info<< "replaced row<1> = " << t3.row<2>() << nl;
+        Info<< "tensor " << t3 << nl;
+
+        triad tr3(t3);
+        Info<< "triad " << tr3 << " :: T() " << tr3.T() << nl;
+    }
+    Info<< nl;
+
 
     tensor t4(3,-2,1,-2,2,0,1, 0, 4);
 
diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
index ce2ee164f2..c522d8f48a 100644
--- a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
+++ b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
@@ -413,7 +413,7 @@ int main(int argc, char *argv[])
     const SubList<face> outsideFaces
     (
         mesh.faces(),
-        mesh.nFaces() - mesh.nInternalFaces(),
+        mesh.nBoundaryFaces(),
         mesh.nInternalFaces()
     );
 
diff --git a/applications/utilities/mesh/conversion/fluentMeshToFoam/fluentMeshToFoam.L b/applications/utilities/mesh/conversion/fluentMeshToFoam/fluentMeshToFoam.L
index 064395df05..e029457d42 100644
--- a/applications/utilities/mesh/conversion/fluentMeshToFoam/fluentMeshToFoam.L
+++ b/applications/utilities/mesh/conversion/fluentMeshToFoam/fluentMeshToFoam.L
@@ -1386,7 +1386,7 @@ int main(int argc, char *argv[])
     DynamicList<polyPatch*> newPatches(nPatches);
 
     // For every boundary face the old patch.
-    labelList facePatchID(pShapeMesh.nFaces()-pShapeMesh.nInternalFaces(), -1);
+    labelList facePatchID(pShapeMesh.nBoundaryFaces(), -1);
     label cMeshFace = pShapeMesh.nInternalFaces();
     label nBoundaries = 0;
 
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index ce5f0bc191..3ab96dd6b0 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -468,7 +468,7 @@ void Foam::backgroundMeshDecomposition::printMeshData
 
     // // Print stats
 
-    // globalIndex globalBoundaryFaces(mesh.nFaces()-mesh.nInternalFaces());
+    // globalIndex globalBoundaryFaces(mesh.nBoundaryFaces());
 
     for (label proci = 0; proci < Pstream::nProcs(); proci++)
     {
@@ -671,7 +671,7 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
         SubList<face>
         (
             mesh_.faces(),
-            mesh_.nFaces() - mesh_.nInternalFaces(),
+            mesh_.nBoundaryFaces(),
             mesh_.nInternalFaces()
         ),
         mesh_.points()
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C
index 2a37049388..7f2ecf8c81 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C
@@ -39,7 +39,7 @@ void Foam::conformalVoronoiMesh::calcNeighbourCellCentres
     pointField& neiCc
 ) const
 {
-    label nBoundaryFaces = mesh.nFaces() - mesh.nInternalFaces();
+    const label nBoundaryFaces = mesh.nBoundaryFaces();
 
     if (neiCc.size() != nBoundaryFaces)
     {
@@ -319,7 +319,7 @@ void Foam::conformalVoronoiMesh::calcFaceZones
     const labelList& faceOwner = mesh.faceOwner();
     const labelList& faceNeighbour = mesh.faceNeighbour();
 
-    labelList neiFaceOwner(mesh.nFaces() - mesh.nInternalFaces(), label(-1));
+    labelList neiFaceOwner(mesh.nBoundaryFaces(), label(-1));
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
@@ -421,7 +421,7 @@ void Foam::conformalVoronoiMesh::calcFaceZones
         )
     );
 
-    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
+    pointField neiCc(mesh.nBoundaryFaces());
     calcNeighbourCellCentres
     (
         mesh,
@@ -498,7 +498,7 @@ void Foam::conformalVoronoiMesh::calcFaceZones
     }
 
 
-//    labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
+//    labelList neiCellSurface(mesh.nBoundaryFaces());
 //
 //    forAll(patches, patchi)
 //    {
diff --git a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
index 578731132f..0364c17774 100644
--- a/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
+++ b/applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C
@@ -119,7 +119,7 @@ void printMeshData(const polyMesh& mesh)
 
     // Print stats
 
-    globalIndex globalBoundaryFaces(mesh.nFaces()-mesh.nInternalFaces());
+    globalIndex globalBoundaryFaces(mesh.nBoundaryFaces());
 
     label maxProcCells = 0;
     label totProcFaces = 0;
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index d899194c19..5daa5766fe 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -525,7 +525,7 @@ void extractSurface
         fileName globalCasePath
         (
             runTime.processorCase()
-          ? runTime.path()/".."/outFileName
+          ? runTime.globalPath()/outFileName
           : runTime.path()/outFileName
         );
         globalCasePath.clean();
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTools.C b/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
index d3d7853436..e230e03990 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTools.C
@@ -314,7 +314,7 @@ void Foam::mergeAndWrite
     }
 
 
-    boolList bndInSet(mesh.nFaces()-mesh.nInternalFaces());
+    boolList bndInSet(mesh.nBoundaryFaces());
     forAll(pbm, patchi)
     {
         const polyPatch& pp = pbm[patchi];
diff --git a/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C b/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C
index 5422183a3f..d6c7160f5d 100644
--- a/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C
+++ b/applications/utilities/mesh/manipulation/mergeOrSplitBaffles/mergeOrSplitBaffles.C
@@ -249,7 +249,7 @@ labelList findBaffles(const polyMesh& mesh, const labelList& boundaryFaces)
         (
             mesh,
             "duplicateFaces",
-            (mesh.nFaces() - mesh.nInternalFaces())/256
+            mesh.nBoundaryFaces()/256
         );
 
         forAll(duplicates, bFacei)
diff --git a/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C b/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C
index 2a05649e77..972503c9ec 100644
--- a/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C
+++ b/applications/utilities/mesh/manipulation/orientFaceZone/orientFaceZone.C
@@ -295,11 +295,7 @@ int main(int argc, char *argv[])
     {
         const polyBoundaryMesh& bm = mesh.boundaryMesh();
 
-        labelList neiStatus
-        (
-            mesh.nFaces()-mesh.nInternalFaces(),
-            orientedSurface::UNVISITED
-        );
+        labelList neiStatus(mesh.nBoundaryFaces(), orientedSurface::UNVISITED);
 
         forAll(faceLabels, i)
         {
diff --git a/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C b/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C
index 88f949782d..a1cbb287f6 100644
--- a/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C
+++ b/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C
@@ -140,7 +140,7 @@ void simpleMarkFeatures
         SubList<face>
         (
             mesh.faces(),
-            mesh.nFaces()-mesh.nInternalFaces(),
+            mesh.nBoundaryFaces(),
             mesh.nInternalFaces()
         ),
         mesh.points()
@@ -238,7 +238,7 @@ void simpleMarkFeatures
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Face centres that need inclusion in the dual mesh
-    labelHashSet featureFaceSet(mesh.nFaces()-mesh.nInternalFaces());
+    labelHashSet featureFaceSet(mesh.nBoundaryFaces());
     // A. boundary faces.
     for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
     {
diff --git a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
index fe0a60070f..2b9ff9c3b8 100644
--- a/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
+++ b/applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C
@@ -353,7 +353,7 @@ void getInterfaceSizes
     // ~~~~~~~~~~~~~~
 
     // Neighbour cellRegion.
-    labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
+    labelList coupledRegion(mesh.nBoundaryFaces());
 
     forAll(coupledRegion, i)
     {
@@ -641,7 +641,7 @@ autoPtr<mapPolyMesh> createRegionMesh
 
 
     // Neighbour cellRegion.
-    labelList coupledRegion(mesh.nFaces()-mesh.nInternalFaces());
+    labelList coupledRegion(mesh.nBoundaryFaces());
 
     forAll(coupledRegion, i)
     {
@@ -1232,7 +1232,7 @@ void getZoneID
     }
 
     // Neighbour zoneID.
-    neiZoneID.setSize(mesh.nFaces()-mesh.nInternalFaces());
+    neiZoneID.setSize(mesh.nBoundaryFaces());
 
     forAll(neiZoneID, i)
     {
@@ -1263,7 +1263,7 @@ void matchRegions
 
     // Get current per cell zoneID
     labelList zoneID(mesh.nCells(), -1);
-    labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+    labelList neiZoneID(mesh.nBoundaryFaces());
     getZoneID(mesh, cellZones, zoneID, neiZoneID);
 
     // Sizes per cellzone
@@ -1552,7 +1552,7 @@ int main(int argc, char *argv[])
     // Existing zoneID
     labelList zoneID(mesh.nCells(), -1);
     // Neighbour zoneID.
-    labelList neiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+    labelList neiZoneID(mesh.nBoundaryFaces());
     getZoneID(mesh, cellZones, zoneID, neiZoneID);
 
 
@@ -1622,7 +1622,7 @@ int main(int argc, char *argv[])
         );
 
         labelList newZoneID(mesh.nCells(), -1);
-        labelList newNeiZoneID(mesh.nFaces()-mesh.nInternalFaces());
+        labelList newNeiZoneID(mesh.nBoundaryFaces());
         getZoneID(mesh, newCellZones, newZoneID, newNeiZoneID);
 
         label unzonedCelli = newZoneID.find(-1);
diff --git a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
index 133161fcc7..664ba34078 100644
--- a/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
+++ b/applications/utilities/parallelProcessing/redistributePar/redistributePar.C
@@ -172,7 +172,7 @@ void printMeshData(const polyMesh& mesh)
 
     // Print stats
 
-    globalIndex globalBoundaryFaces(mesh.nFaces()-mesh.nInternalFaces());
+    globalIndex globalBoundaryFaces(mesh.nBoundaryFaces());
 
     label maxProcCells = 0;
     label totProcFaces = 0;
@@ -2451,7 +2451,7 @@ int main(int argc, char *argv[])
             << " on all processors" << nl << endl;
         instantList timeDirs;
 
-        const fileName basePath(args.rootPath()/args.globalCaseName());
+        const fileName basePath(args.globalPath());
 
         if (Pstream::master())
         {
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
index 3de907152b..8161e4ae15 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsight/foamToEnsight.C
@@ -225,7 +225,7 @@ int main(int argc, char *argv[])
     fileName ensightDir = args.lookupOrDefault<word>("name", "EnSight");
     if (!ensightDir.isAbsolute())
     {
-        ensightDir = args.rootPath()/args.globalCaseName()/ensightDir;
+        ensightDir = args.globalPath()/ensightDir;
     }
 
 
diff --git a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C
index 13adf3b9fd..01c9d0244b 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToEnsightParts/foamToEnsightParts.C
@@ -194,7 +194,7 @@ int main(int argc, char *argv[])
     fileName ensightDir = args.lookupOrDefault<word>("name", "Ensight");
     if (!ensightDir.isAbsolute())
     {
-        ensightDir = args.rootPath()/args.globalCaseName()/ensightDir;
+        ensightDir = args.globalPath()/ensightDir;
     }
 
     //
diff --git a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
index b90e4cefea..ddefd7f7ee 100644
--- a/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
+++ b/applications/utilities/postProcessing/miscellaneous/postChannel/channelIndex.C
@@ -55,7 +55,7 @@ void Foam::channelIndex::walkOppositeFaces
 {
     const cellList& cells = mesh.cells();
     const faceList& faces = mesh.faces();
-    label nBnd = mesh.nFaces() - mesh.nInternalFaces();
+    const label nBnd = mesh.nBoundaryFaces();
 
     DynamicList<label> frontFaces(startFaces);
     forAll(frontFaces, i)
diff --git a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
index e48f7b48e3..ff77af6bb6 100644
--- a/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
+++ b/applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C
@@ -128,7 +128,7 @@ int main(int argc, char *argv[])
     }
 
     // Sync agglomeration across coupled patches
-    labelList nbrAgglom(mesh.nFaces() - mesh.nInternalFaces(), -1);
+    labelList nbrAgglom(mesh.nBoundaryFaces(), -1);
 
     forAll(boundary, patchi)
     {
diff --git a/applications/utilities/preProcessing/setFields/setFields.C b/applications/utilities/preProcessing/setFields/setFields.C
index f7578ee9de..554f8e1ac9 100644
--- a/applications/utilities/preProcessing/setFields/setFields.C
+++ b/applications/utilities/preProcessing/setFields/setFields.C
@@ -239,7 +239,7 @@ bool setFaceFieldType
         const Type& value = pTraits<Type>(fieldValueStream);
 
         // Create flat list of selected faces and their value.
-        Field<Type> allBoundaryValues(mesh.nFaces()-mesh.nInternalFaces());
+        Field<Type> allBoundaryValues(mesh.nBoundaryFaces());
         forAll(field.boundaryField(), patchi)
         {
             SubField<Type>
@@ -447,7 +447,7 @@ int main(int argc, char *argv[])
             (
                 mesh,
                 "faceSet",
-                (mesh.nFaces()-mesh.nInternalFaces())/10+1
+                mesh.nBoundaryFaces()/10+1
             );
 
             source->applyToSet
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 7d575efa82..3670eceab1 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -82,10 +82,7 @@ triSurface triangulate
     const polyMesh& mesh = bMesh.mesh();
 
     // Storage for surfaceMesh. Size estimate.
-    DynamicList<labelledTri> triangles
-    (
-        mesh.nFaces() - mesh.nInternalFaces()
-    );
+    DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
 
     label newPatchI = 0;
     label localTriFaceI = 0;
diff --git a/applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C b/applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C
index f18cf71f50..a3eedaa734 100644
--- a/applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C
+++ b/applications/utilities/surface/surfaceMeshTriangulate/surfaceMeshTriangulate.C
@@ -411,7 +411,7 @@ int main(int argc, char *argv[])
               ? outFileName
               : (
                     runTime.processorCase()
-                  ? runTime.rootPath()/runTime.globalCaseName()/outFileName
+                  ? runTime.globalPath()/outFileName
                   : runTime.path()/outFileName
                 )
             );
diff --git a/src/OpenFOAM/db/Time/TimePaths.H b/src/OpenFOAM/db/Time/TimePaths.H
index 38189a4892..c5067b9227 100644
--- a/src/OpenFOAM/db/Time/TimePaths.H
+++ b/src/OpenFOAM/db/Time/TimePaths.H
@@ -64,7 +64,7 @@ class TimePaths
         const word constant_;
 
 
-    // Private Member functions
+    // Private Member Functions
 
         //- Determine from case name whether it is a processor directory
         bool detectProcessorCase();
@@ -105,7 +105,7 @@ public:
         );
 
 
-    // Member functions
+    // Member Functions
 
         //- Return true if this is a processor case
         bool processorCase() const
@@ -164,12 +164,18 @@ public:
         //- \c ../constant() for parallel runs.
         fileName caseConstant() const;
 
-        //- Return path
+        //- Return path for the case
         fileName path() const
         {
             return rootPath()/caseName();
         }
 
+        //- Return global path for the case
+        fileName globalPath() const
+        {
+            return rootPath()/globalCaseName();
+        }
+
         //- Return system path
         fileName systemPath() const
         {
diff --git a/src/OpenFOAM/fields/Fields/transformList/transformList.C b/src/OpenFOAM/fields/Fields/transformList/transformList.C
index 87f8c05cc1..7088bcd7e0 100644
--- a/src/OpenFOAM/fields/Fields/transformList/transformList.C
+++ b/src/OpenFOAM/fields/Fields/transformList/transformList.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,14 +34,14 @@ Foam::List<T> Foam::transform
     const UList<T>& field
 )
 {
-    List<T> newField(field.size());
+    List<T> result(field.size());
 
     forAll(field, i)
     {
-        newField[i] = transform(rotTensor, field[i]);
+        result[i] = transform(rotTensor, field[i]);
     }
 
-    return newField;
+    return result;
 }
 
 
@@ -60,10 +60,7 @@ void Foam::transformList(const tensorField& rotTensor, UList<T>& field)
 {
     if (rotTensor.size() == 1)
     {
-        forAll(field, i)
-        {
-            field[i] = transform(rotTensor[0], field[i]);
-        }
+        transformList(rotTensor[0], field);
     }
     else if (rotTensor.size() == field.size())
     {
@@ -85,9 +82,10 @@ void Foam::transformList(const tensorField& rotTensor, UList<T>& field)
 template<class T>
 void Foam::transformList(const tensor& rotTensor, Map<T>& field)
 {
-    forAllIter(typename Map<T>, field, iter)
+    forAllIters(field, iter)
     {
-        iter() = transform(rotTensor[0], iter());
+        T& value = iter.object();
+        value = transform(rotTensor, value);
     }
 }
 
@@ -97,10 +95,7 @@ void Foam::transformList(const tensorField& rotTensor, Map<T>& field)
 {
     if (rotTensor.size() == 1)
     {
-        forAllIter(typename Map<T>, field, iter)
-        {
-            iter() = transform(rotTensor[0], iter());
-        }
+        transformList(rotTensor[0], field);
     }
     else
     {
@@ -115,9 +110,10 @@ void Foam::transformList(const tensorField& rotTensor, Map<T>& field)
 template<class T>
 void Foam::transformList(const tensor& rotTensor, EdgeMap<T>& field)
 {
-    forAllIter(typename EdgeMap<T>, field, iter)
+    forAllIters(field, iter)
     {
-        iter() = transform(rotTensor[0], iter());
+        T& value = iter.object();
+        value = transform(rotTensor, value);
     }
 }
 
@@ -127,10 +123,7 @@ void Foam::transformList(const tensorField& rotTensor, EdgeMap<T>& field)
 {
     if (rotTensor.size() == 1)
     {
-        forAllIter(typename EdgeMap<T>, field, iter)
-        {
-            iter() = transform(rotTensor[0], iter());
-        }
+        transformList(rotTensor[0], field);
     }
     else
     {
diff --git a/src/OpenFOAM/fields/Fields/transformList/transformList.H b/src/OpenFOAM/fields/Fields/transformList/transformList.H
index 8b485bec41..9f827d0928 100644
--- a/src/OpenFOAM/fields/Fields/transformList/transformList.H
+++ b/src/OpenFOAM/fields/Fields/transformList/transformList.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,7 +25,7 @@ InClass
     Foam
 
 Description
-    Spatial transformation functions for primitive fields.
+    Spatial transformation functions for list of values and primitive fields.
 
 SourceFiles
     transformList.C
@@ -41,7 +41,6 @@ SourceFiles
 #include "edgeHashes.H"
 #include "tensorField.H"
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -49,93 +48,63 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//- Extend transform to work on list.
+//- Apply transform to a list of elements, returning a copy.
 template<class T>
-List<T> transform
-(
-    const tensor& rotTensor,
-    const UList<T>& field
-);
+List<T> transform(const tensor& rotTensor, const UList<T>& field);
 
-//- Apply transformation to list. Either single transformation tensor
-//  or one tensor per element.
+//- Inplace transform a list of elements.
 template<class T>
-void transformList(const tensor&, UList<T>&);
+void transformList(const tensor& rotTensor, UList<T>& field);
+
+//- Inplace transform a list of elements using one tensor per element.
 template<class T>
-void transformList(const tensorField&, UList<T>&);
+void transformList(const tensorField& rotTensor, UList<T>& field);
 
+//- Inplace transform a Map of elements.
 template<class T>
-void transformList(const tensor&, Map<T>&);
+void transformList(const tensor& rotTensor, Map<T>& field);
+
+//- Inplace transform a Map of elements using one tensor per element.
+//  Using multiple tensors is ill-defined (Fatal).
 template<class T>
-void transformList(const tensorField&, Map<T>&);
+void transformList(const tensorField& rotTensor, Map<T>& field);
 
+//- Inplace transform a Map of elements.
+//  Using multiple tensors is ill-defined (Fatal).
 template<class T>
-void transformList(const tensor&, EdgeMap<T>&);
+void transformList(const tensor& rotTensor, EdgeMap<T>& field);
+
+//- Inplace transform a Map of elements using one tensor per element.
+//  Using multiple tensors is ill-defined (Fatal).
 template<class T>
-void transformList(const tensorField&, EdgeMap<T>&);
+void transformList(const tensorField& rotTensor, EdgeMap<T>& field);
 
 
-// Specialisations for bool
-template<>
-inline void transformList(const tensor&, UList<bool>&)
-{}
-template<>
-inline void transformList(const tensorField&, UList<bool>&)
-{}
-template<>
-inline void transformList(const tensor&, Map<bool>&)
-{}
-template<>
-inline void transformList(const tensorField&, Map<bool>&)
-{}
-template<>
-inline void transformList(const tensor&, EdgeMap<bool>&)
-{}
-template<>
-inline void transformList(const tensorField&, EdgeMap<bool>&)
-{}
+// Specialisations for bool (no-op)
+template<> inline void transformList(const tensor&, UList<bool>&) {}
+template<> inline void transformList(const tensorField&, UList<bool>&) {}
+template<> inline void transformList(const tensor&, Map<bool>&) {}
+template<> inline void transformList(const tensorField&, Map<bool>&) {}
+template<> inline void transformList(const tensor&, EdgeMap<bool>&) {}
+template<> inline void transformList(const tensorField&, EdgeMap<bool>&) {}
 
 
-// Specialisations for label
-template<>
-inline void transformList(const tensor&, labelUList&)
-{}
-template<>
-inline void transformList(const tensorField&, labelUList&)
-{}
-template<>
-inline void transformList(const tensor&, Map<label>&)
-{}
-template<>
-inline void transformList(const tensorField&, Map<label>&)
-{}
-template<>
-inline void transformList(const tensor&, EdgeMap<label>&)
-{}
-template<>
-inline void transformList(const tensorField&, EdgeMap<label>&)
-{}
+// Specialisations for label (no-op)
+template<> inline void transformList(const tensor&, labelUList&) {}
+template<> inline void transformList(const tensorField&, labelUList&) {}
+template<> inline void transformList(const tensor&, Map<label>&) {}
+template<> inline void transformList(const tensorField&, Map<label>&) {}
+template<> inline void transformList(const tensor&, EdgeMap<label>&) {}
+template<> inline void transformList(const tensorField&, EdgeMap<label>&) {}
 
 
-// Specialisations for scalar
-template<>
-inline void transformList(const tensor&, UList<scalar>&)
-{}
-template<>
-inline void transformList(const tensorField&, UList<scalar>&)
-{}
-template<>
-inline void transformList(const tensor&, Map<scalar>&)
-{}
-template<>
-inline void transformList(const tensorField&, Map<scalar>&)
-{}
-template<>
-inline void transformList(const tensor&, EdgeMap<scalar>&)
-{}
-template<>
-inline void transformList(const tensorField&, EdgeMap<scalar>&)
-{}
+// Specialisations for scalar (no-op)
+template<> inline void transformList(const tensor&, UList<scalar>&) {}
+template<> inline void transformList(const tensorField&, UList<scalar>&) {}
+template<> inline void transformList(const tensor&, Map<scalar>&) {}
+template<> inline void transformList(const tensorField&, Map<scalar>&) {}
+template<> inline void transformList(const tensor&, EdgeMap<scalar>&) {}
+template<> inline void transformList(const tensorField&, EdgeMap<scalar>&) {}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/global/argList/argList.H b/src/OpenFOAM/global/argList/argList.H
index 42448c3b5c..ba1baa8f3a 100644
--- a/src/OpenFOAM/global/argList/argList.H
+++ b/src/OpenFOAM/global/argList/argList.H
@@ -293,7 +293,7 @@ public:
         //- Return case name (parallel run) or global case (serial run)
         inline const fileName& caseName() const;
 
-        //- Return case name
+        //- Return global case name
         inline const fileName& globalCaseName() const;
 
         //- Return parRunControl
@@ -302,6 +302,9 @@ public:
         //- Return the path to the caseName
         inline fileName path() const;
 
+        //- Return the path to the globalCaseName
+        inline fileName globalPath() const;
+
         //- Return the number of arguments
         inline label size() const;
 
diff --git a/src/OpenFOAM/global/argList/argListI.H b/src/OpenFOAM/global/argList/argListI.H
index 83b521f923..8457a4cf82 100644
--- a/src/OpenFOAM/global/argList/argListI.H
+++ b/src/OpenFOAM/global/argList/argListI.H
@@ -87,6 +87,12 @@ inline Foam::fileName Foam::argList::path() const
 }
 
 
+inline Foam::fileName Foam::argList::globalPath() const
+{
+    return rootPath()/globalCaseName();
+}
+
+
 inline Foam::label Foam::argList::size() const
 {
     return args_.size();
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index ae6d7f28a2..feb7781a19 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -1325,7 +1325,7 @@ void Foam::globalMeshData::calcGlobalPointBoundaryFaces() const
     // Global indices for boundary faces
     globalBoundaryFaceNumberingPtr_.reset
     (
-        new globalIndex(mesh_.nFaces()-mesh_.nInternalFaces())
+        new globalIndex(mesh_.nBoundaryFaces())
     );
     globalIndex& globalIndices = globalBoundaryFaceNumberingPtr_();
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
index 9182fb5a90..741aeef3e8 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C
@@ -450,10 +450,7 @@ const Foam::labelList& Foam::polyBoundaryMesh::patchID() const
 {
     if (!patchIDPtr_.valid())
     {
-        patchIDPtr_.reset
-        (
-            new labelList(mesh_.nFaces() - mesh_.nInternalFaces())
-        );
+        patchIDPtr_.reset(new labelList(mesh_.nBoundaryFaces()));
         labelList& list = *patchIDPtr_;
 
         const polyPatchList& patches = *this;
@@ -613,13 +610,21 @@ Foam::wordList Foam::polyBoundaryMesh::physicalTypes() const
 }
 
 
+Foam::label Foam::polyBoundaryMesh::start() const
+{
+    return mesh_.nInternalFaces();
+}
+
+
+Foam::label Foam::polyBoundaryMesh::nFaces() const
+{
+    return mesh_.nBoundaryFaces();
+}
+
+
 Foam::labelRange Foam::polyBoundaryMesh::range() const
 {
-    return labelRange
-    (
-        mesh_.nInternalFaces(),
-        mesh_.nFaces() - mesh_.nInternalFaces()
-    );
+    return labelRange(mesh_.nInternalFaces(), mesh_.nBoundaryFaces());
 }
 
 
diff --git a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
index 0947df807d..c97da3db7b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.H
@@ -153,7 +153,7 @@ public:
         //  polyBoundaryMesh and not parallel
         const List<labelPairList>& neighbourEdges() const;
 
-        //- Return the number of non-processor patches
+        //- The number of patches before the first processor patch.
         label nNonProcessor() const;
 
         //- Return a list of patch names
@@ -165,8 +165,16 @@ public:
         //- Return a list of physical types
         wordList physicalTypes() const;
 
-        //- Return the range used for all boundary faces
-        //  Spans [nInternalFaces, nFaces)
+        //- The start label of the boundary faces in the polyMesh face list
+        //  Same as mesh.nInternalFaces()
+        label start() const;
+
+        //- The number of boundary faces in the underlying mesh
+        //  Same as mesh.nBoundaryFaces()
+        label nFaces() const;
+
+        //- The face range for all boundary faces
+        //  Spans [nInternalFaces, nFaces) of the underlying mesh
         labelRange range() const;
 
         //- Return the range used for boundary faces on patchi.
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
index f2389857d9..a84f56d5d7 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.C
@@ -377,7 +377,7 @@ bool Foam::polyMeshTetDecomposition::checkFaceTets
     const vectorField& fc = mesh.faceCentres();
 
     // Calculate coupled cell centre
-    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
+    pointField neiCc(mesh.nBoundaryFaces());
 
     for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei)
     {
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
index d48d7e30ca..f1cb956a57 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -293,6 +293,12 @@ Foam::wordList Foam::polyPatch::constraintTypes()
 }
 
 
+Foam::label Foam::polyPatch::offset() const
+{
+    return start_ - boundaryMesh().start();
+}
+
+
 const Foam::polyBoundaryMesh& Foam::polyPatch::boundaryMesh() const
 {
     return boundaryMesh_;
@@ -330,6 +336,25 @@ Foam::tmp<Foam::vectorField> Foam::polyPatch::faceCellCentres() const
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::polyPatch::areaFraction() const
+{
+    tmp<scalarField> tfraction(new scalarField(size()));
+    scalarField& fraction = tfraction.ref();
+
+    const vectorField::subField faceAreas = this->faceAreas();
+    const pointField& points = this->points();
+
+    forAll(*this, facei)
+    {
+        const face& curFace = this->operator[](facei);
+        fraction[facei] =
+            mag(faceAreas[facei])/(curFace.mag(points) + ROOTVSMALL);
+    }
+
+    return tfraction;
+}
+
+
 const Foam::labelUList& Foam::polyPatch::faceCells() const
 {
     if (!faceCellsPtr_)
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
index a43c9e9844..d661304d52 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,10 +49,7 @@ SourceFiles
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Forward declaration of friend functions and operators
-
+// Forward declarations
 class polyBoundaryMesh;
 class polyPatch;
 class PstreamBuffers;
@@ -304,6 +301,10 @@ public:
 
     // Member Functions
 
+        //- The offset where this patch starts in the boundary face list
+        //  The value is the same as  patch.start() - mesh.nInternalFaces()
+        label offset() const;
+
         //- Return start label of this patch in the polyMesh face list
         label start() const
         {
@@ -372,6 +373,10 @@ public:
             //- Return face cell centres
             tmp<vectorField> faceCellCentres() const;
 
+            //- Return the area fraction as the ratio of the stored face area
+            //- and the area given by the face points
+            tmp<scalarField> areaFraction() const;
+
 
         // Addressing into mesh
 
diff --git a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.C b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.C
index 99bf41e4c3..1437834a46 100644
--- a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.C
+++ b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.C
@@ -44,7 +44,7 @@ void Foam::syncTools::swapBoundaryCellPositions
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    neighbourCellData.resize(mesh.nFaces()-mesh.nInternalFaces());
+    neighbourCellData.resize(mesh.nBoundaryFaces());
 
     for (const polyPatch& pp : patches)
     {
diff --git a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
index 68f47b1494..7f32972586 100644
--- a/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
+++ b/src/OpenFOAM/meshes/polyMesh/syncTools/syncTools.H
@@ -395,7 +395,7 @@ public:
                 SubList<T> bndValues
                 (
                     faceValues,
-                    mesh.nFaces()-mesh.nInternalFaces(),
+                    mesh.nBoundaryFaces(),
                     mesh.nInternalFaces()
                 );
 
@@ -420,7 +420,7 @@ public:
                 SubList<point> bndValues
                 (
                     positions,
-                    mesh.nFaces()-mesh.nInternalFaces(),
+                    mesh.nBoundaryFaces(),
                     mesh.nInternalFaces()
                 );
                 syncBoundaryFaceList
@@ -476,7 +476,7 @@ public:
                 SubList<T> bndValues
                 (
                     faceValues,
-                    mesh.nFaces()-mesh.nInternalFaces(),
+                    mesh.nBoundaryFaces(),
                     mesh.nInternalFaces()
                 );
                 syncBoundaryFaceList
diff --git a/src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C b/src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C
index e96be7fe10..c0b15c78da 100644
--- a/src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C
+++ b/src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C
@@ -929,7 +929,7 @@ void Foam::syncTools::syncBoundaryFaceList
     const bool parRun
 )
 {
-    const label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
+    const label nBFaces = mesh.nBoundaryFaces();
 
     if (faceValues.size() != nBFaces)
     {
@@ -1160,7 +1160,7 @@ void Foam::syncTools::swapBoundaryCellList
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-    neighbourCellData.resize(mesh.nFaces()-mesh.nInternalFaces());
+    neighbourCellData.resize(mesh.nBoundaryFaces());
 
     for (const polyPatch& pp : patches)
     {
diff --git a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
index 8bf7480557..d7543f4fd5 100644
--- a/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
+++ b/src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C
@@ -537,8 +537,8 @@ bool Foam::faceZone::checkParallelSync(const bool report) const
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     {
-        boolList neiZoneFace(mesh.nFaces()-mesh.nInternalFaces(), false);
-        boolList neiZoneFlip(mesh.nFaces()-mesh.nInternalFaces(), false);
+        boolList neiZoneFace(mesh.nBoundaryFaces(), false);
+        boolList neiZoneFlip(mesh.nBoundaryFaces(), false);
         forAll(*this, i)
         {
             const label facei = operator[](i);
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index 5e1f71275a..cbd7df5f90 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -91,7 +91,7 @@ class primitiveMesh
             //- Number of internal edges using 0 or 1 boundary points
             mutable label nInternal1Edges_;
 
-            //- Number of internal edges using 0,1 or 2boundary points
+            //- Number of internal edges using 0,1 or 2 boundary points
             mutable label nInternalEdges_;
 
             //- Number of edges
@@ -478,24 +478,39 @@ public:
 
             // Mesh size parameters
 
+                //- Number of mesh points
                 inline label nPoints() const;
+
+                //- Number of mesh edges
                 inline label nEdges() const;
-                inline label nInternalFaces() const;
+
+                //- Number of mesh faces
                 inline label nFaces() const;
+
+                //- Number of mesh cells
                 inline label nCells() const;
 
-                // If points are ordered (nInternalPoints != -1):
+                //- Number of internal faces
+                inline label nInternalFaces() const;
 
-                    //- Points not on boundary
-                    inline label nInternalPoints() const;
+                //- Number of boundary faces (== nFaces - nInternalFaces)
+                inline label nBoundaryFaces() const;
 
-                    //- Internal edges (i.e. not on boundary face) using
-                    //  no boundary point
-                    inline label nInternal0Edges() const;
-                    //- Internal edges using 0 or 1 boundary point
-                    inline label nInternal1Edges() const;
-                    //- Internal edges using 0,1 or 2 boundary points
-                    inline label nInternalEdges() const;
+
+            // If points are ordered (nInternalPoints != -1):
+
+                //- Points not on boundary
+                inline label nInternalPoints() const;
+
+                //- Internal edges (i.e. not on boundary face) using
+                //- no boundary point
+                inline label nInternal0Edges() const;
+
+                //- Internal edges using 0 or 1 boundary point
+                inline label nInternal1Edges() const;
+
+                //- Internal edges using 0,1 or 2 boundary points
+                inline label nInternalEdges() const;
 
 
             // Primitive mesh data
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H
index 42ad4c2687..df14e3c362 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -21,30 +21,23 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label primitiveMesh::nInternalPoints() const
+inline Foam::label Foam::primitiveMesh::nInternalPoints() const
 {
     return nInternalPoints_;
 }
 
 
-inline label primitiveMesh::nPoints() const
+inline Foam::label Foam::primitiveMesh::nPoints() const
 {
     return nPoints_;
 }
 
 
-inline label primitiveMesh::nInternal0Edges() const
+inline Foam::label Foam::primitiveMesh::nInternal0Edges() const
 {
     // Force edge calculation
     (void)nEdges();
@@ -52,7 +45,7 @@ inline label primitiveMesh::nInternal0Edges() const
 }
 
 
-inline label primitiveMesh::nInternal1Edges() const
+inline Foam::label Foam::primitiveMesh::nInternal1Edges() const
 {
     // Force edge calculation
     (void)nEdges();
@@ -60,7 +53,7 @@ inline label primitiveMesh::nInternal1Edges() const
 }
 
 
-inline label primitiveMesh::nInternalEdges() const
+inline Foam::label Foam::primitiveMesh::nInternalEdges() const
 {
     // Force edge calculation
     (void)nEdges();
@@ -68,7 +61,7 @@ inline label primitiveMesh::nInternalEdges() const
 }
 
 
-inline label primitiveMesh::nEdges() const
+inline Foam::label Foam::primitiveMesh::nEdges() const
 {
     if (nEdges_ < 0)
     {
@@ -79,134 +72,136 @@ inline label primitiveMesh::nEdges() const
 }
 
 
-inline label primitiveMesh::nInternalFaces() const
+inline Foam::label Foam::primitiveMesh::nInternalFaces() const
 {
     return nInternalFaces_;
 }
 
 
-inline label primitiveMesh::nFaces() const
+inline Foam::label Foam::primitiveMesh::nBoundaryFaces() const
+{
+    return (nFaces_ - nInternalFaces_);
+}
+
+
+inline Foam::label Foam::primitiveMesh::nFaces() const
 {
     return nFaces_;
 }
 
 
-inline label primitiveMesh::nCells() const
+inline Foam::label Foam::primitiveMesh::nCells() const
 {
     return nCells_;
 }
 
 
-inline bool primitiveMesh::isInternalFace(const label faceIndex) const
+inline bool Foam::primitiveMesh::isInternalFace(const label faceIndex) const
 {
-    return faceIndex < nInternalFaces();
+    return faceIndex < nInternalFaces_;
 }
 
 
-inline bool primitiveMesh::hasCellShapes() const
+inline bool Foam::primitiveMesh::hasCellShapes() const
 {
     return cellShapesPtr_;
 }
 
 
-inline bool primitiveMesh::hasEdges() const
+inline bool Foam::primitiveMesh::hasEdges() const
 {
     return edgesPtr_;
 }
 
 
-inline bool primitiveMesh::hasCellCells() const
+inline bool Foam::primitiveMesh::hasCellCells() const
 {
     return ccPtr_;
 }
 
 
-inline bool primitiveMesh::hasEdgeCells() const
+inline bool Foam::primitiveMesh::hasEdgeCells() const
 {
     return ecPtr_;
 }
 
 
-inline bool primitiveMesh::hasPointCells() const
+inline bool Foam::primitiveMesh::hasPointCells() const
 {
     return pcPtr_;
 }
 
 
-inline bool primitiveMesh::hasCells() const
+inline bool Foam::primitiveMesh::hasCells() const
 {
     return cfPtr_;
 }
 
 
-inline bool primitiveMesh::hasEdgeFaces() const
+inline bool Foam::primitiveMesh::hasEdgeFaces() const
 {
     return efPtr_;
 }
 
 
-inline bool primitiveMesh::hasPointFaces() const
+inline bool Foam::primitiveMesh::hasPointFaces() const
 {
     return pfPtr_;
 }
 
 
-inline bool primitiveMesh::hasCellEdges() const
+inline bool Foam::primitiveMesh::hasCellEdges() const
 {
     return cePtr_;
 }
 
 
-inline bool primitiveMesh::hasFaceEdges() const
+inline bool Foam::primitiveMesh::hasFaceEdges() const
 {
     return fePtr_;
 }
 
 
-inline bool primitiveMesh::hasPointEdges() const
+inline bool Foam::primitiveMesh::hasPointEdges() const
 {
     return pePtr_;
 }
 
 
-inline bool primitiveMesh::hasPointPoints() const
+inline bool Foam::primitiveMesh::hasPointPoints() const
 {
     return ppPtr_;
 }
 
 
-inline bool primitiveMesh::hasCellPoints() const
+inline bool Foam::primitiveMesh::hasCellPoints() const
 {
     return cpPtr_;
 }
 
 
-inline bool primitiveMesh::hasCellCentres() const
+inline bool Foam::primitiveMesh::hasCellCentres() const
 {
     return cellCentresPtr_;
 }
 
 
-inline bool primitiveMesh::hasFaceCentres() const
+inline bool Foam::primitiveMesh::hasFaceCentres() const
 {
     return faceCentresPtr_;
 }
 
 
-inline bool primitiveMesh::hasCellVolumes() const
+inline bool Foam::primitiveMesh::hasCellVolumes() const
 {
     return cellVolumesPtr_;
 }
 
 
-inline bool primitiveMesh::hasFaceAreas() const
+inline bool Foam::primitiveMesh::hasFaceAreas() const
 {
     return faceAreasPtr_;
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
 // ************************************************************************* //
diff --git a/src/OpenFOAM/primitives/Tensor/Tensor.H b/src/OpenFOAM/primitives/Tensor/Tensor.H
index ee670f07b9..d39f8b10e1 100644
--- a/src/OpenFOAM/primitives/Tensor/Tensor.H
+++ b/src/OpenFOAM/primitives/Tensor/Tensor.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -95,22 +95,22 @@ public:
 
         //- Construct given MatrixSpace of the same rank
         template<class Cmpt2>
-        inline Tensor(const MatrixSpace<Tensor<Cmpt2>, Cmpt2, 3, 3>&);
+        inline Tensor(const MatrixSpace<Tensor<Cmpt2>, Cmpt2, 3, 3>& vs);
 
         //- Construct given VectorSpace of the same rank
         template<class Cmpt2>
-        inline Tensor(const VectorSpace<Tensor<Cmpt2>, Cmpt2, 9>&);
+        inline Tensor(const VectorSpace<Tensor<Cmpt2>, Cmpt2, 9>& vs);
 
         //- Construct given SphericalTensor
-        inline Tensor(const SphericalTensor<Cmpt>&);
+        inline Tensor(const SphericalTensor<Cmpt>& st);
 
         //- Construct given SymmTensor
-        inline Tensor(const SymmTensor<Cmpt>&);
+        inline Tensor(const SymmTensor<Cmpt>& st);
 
-        //- Construct given triad
-        inline Tensor(const Vector<Vector<Cmpt>>&);
+        //- Construct given triad of row vectors
+        inline Tensor(const Vector<Vector<Cmpt>>& tr);
 
-        //- Construct given the three vector components
+        //- Construct given the three row vectors
         inline Tensor
         (
             const Vector<Cmpt>& x,
@@ -139,7 +139,7 @@ public:
         );
 
         //- Construct from Istream
-        inline Tensor(Istream&);
+        inline Tensor(Istream& is);
 
 
     // Member Functions
@@ -166,12 +166,72 @@ public:
             inline Cmpt& zy();
             inline Cmpt& zz();
 
+
+        // Column-vector access.
+
+            //- Extract vector for column 0
+            inline Vector<Cmpt> cx() const;
+
+            //- Extract vector for column 1
+            inline Vector<Cmpt> cy() const;
+
+            //- Extract vector for column 2
+            inline Vector<Cmpt> cz() const;
+
+            //- Extract vector for given column.
+            //  Compile-time check of column index.
+            template<direction Col>
+            inline Vector<Cmpt> col() const;
+
+            //- Extract vector for given column (0,1,2).
+            //  Runtime check of column index.
+            inline Vector<Cmpt> col(const direction c) const;
+
+            //- Set values of given column
+            //  Compile-time check of column index.
+            template<direction Col>
+            inline void col(const Vector<Cmpt>& v);
+
+            //- Set values of given column (0,1,2)
+            //  Runtime check of column index.
+            inline void col(const direction c, const Vector<Cmpt>& v);
+
+
         // Row-vector access.
 
+            //- Extract vector for row 0
             inline Vector<Cmpt> x() const;
+
+            //- Extract vector for row 1
             inline Vector<Cmpt> y() const;
+
+            //- Extract vector for row 2
             inline Vector<Cmpt> z() const;
-            inline Vector<Cmpt> vectorComponent(const direction) const;
+
+            //- Extract vector for given row.
+            //  Compile-time check of row index.
+            template<direction Row>
+            inline Vector<Cmpt> row() const;
+
+            //- Extract vector for given row (0,1,2)
+            //  Runtime check of row index.
+            inline Vector<Cmpt> row(const direction r) const;
+
+            //- Set values of given row
+            //  Compile-time check of row index.
+            template<direction Row>
+            inline void row(const Vector<Cmpt>& v);
+
+            //- Set values of given row (0,1,2)
+            //  Runtime check of row index.
+            inline void row(const direction r, const Vector<Cmpt>& v);
+
+            //- Return vector for given row (0,1,2)
+            //  Runtime check of row index.
+            inline Vector<Cmpt> vectorComponent(const direction cmpt) const;
+
+
+    // Tensor Operations
 
         //- Return transpose
         inline Tensor<Cmpt> T() const;
@@ -179,11 +239,22 @@ public:
         //- Return inverse
         inline Tensor<Cmpt> inv() const;
 
+        //- Inner-product of this with another Tensor.
+        inline Tensor<Cmpt> inner(const Tensor<Cmpt>& t2) const;
+
+        //- Inner-product of this with transpose of another Tensor.
+        //  Primarily useful for coordinate transformations
+        //  (where transpose is the same as the inverse).
+        inline Tensor<Cmpt> innerT(const Tensor<Cmpt>& t2) const;
+
+        //- Schur-product of this with another Tensor.
+        inline Tensor<Cmpt> schur(const Tensor<Cmpt>& t2) const;
+
 
     // Member Operators
 
-        //- Inner-product with a Tensor
-        inline void operator&=(const Tensor<Cmpt>&);
+        //- Assign inner-product of this with another Tensor.
+        inline void operator&=(const Tensor<Cmpt>& t);
 
         //- Inherit MatrixSpace assignment operators
         using Tensor::msType::operator=;
@@ -198,7 +269,7 @@ public:
         //- Assign to a SymmTensor
         inline void operator=(const SymmTensor<Cmpt>&);
 
-        //- Assign to a triad
+        //- Assign to a triad of row vectors
         inline void operator=(const Vector<Vector<Cmpt>>&);
 };
 
diff --git a/src/OpenFOAM/primitives/Tensor/TensorI.H b/src/OpenFOAM/primitives/Tensor/TensorI.H
index c8ecbc7ce6..933b39a5c2 100644
--- a/src/OpenFOAM/primitives/Tensor/TensorI.H
+++ b/src/OpenFOAM/primitives/Tensor/TensorI.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -296,32 +296,179 @@ inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::z() const
 }
 
 
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::cx() const
+{
+    return Vector<Cmpt>(this->v_[XX], this->v_[YX], this->v_[ZX]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::cy() const
+{
+    return Vector<Cmpt>(this->v_[XY], this->v_[YY], this->v_[ZY]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::cz() const
+{
+    return Vector<Cmpt>(this->v_[XZ], this->v_[YZ], this->v_[ZZ]);
+}
+
+
+template<class Cmpt>
+template<Foam::direction Col>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::col() const
+{
+    if (Col == 0) return cx();
+    else if (Col == 1) return cy();
+    else if (Col == 2) return cz();
+
+    static_assert(Col < 3, "Invalid column access");
+    return Zero;
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::col(const direction c) const
+{
+    switch (c)
+    {
+        case 0: return cx(); break;
+        case 1: return cy(); break;
+        case 2: return cz(); break;
+    }
+
+    FatalErrorInFunction
+        << "Invalid column access " << c << abort(FatalError);
+
+    return Zero;
+}
+
+
+template<class Cmpt>
+template<Foam::direction Row>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::row() const
+{
+    if (Row == 0) return x();
+    else if (Row == 1) return y();
+    else if (Row == 2) return z();
+
+    static_assert(Row < 3, "Invalid row access");
+    return Zero;
+}
+
+
+template<class Cmpt>
+inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::row(const direction r) const
+{
+    switch (r)
+    {
+        case 0: return x(); break;
+        case 1: return y(); break;
+        case 2: return z(); break;
+    }
+
+    FatalErrorInFunction
+        << "Invalid row access " << r << abort(FatalError);
+
+    return Zero;
+}
+
+
+template<class Cmpt>
+template<Foam::direction Col>
+inline void Foam::Tensor<Cmpt>::col(const Vector<Cmpt>& v)
+{
+    if (Col == 0)
+    {
+        this->v_[XX] = v.x();
+        this->v_[YX] = v.y();
+        this->v_[ZX] = v.z();
+    }
+    else if (Col == 1)
+    {
+        this->v_[XY] = v.x();
+        this->v_[YY] = v.y();
+        this->v_[ZY] = v.z();
+    }
+    else if (Col == 2)
+    {
+        this->v_[XZ] = v.x();
+        this->v_[YZ] = v.y();
+        this->v_[ZZ] = v.z();
+    }
+
+    static_assert(Col < 3, "Invalid column access");
+}
+
+
+template<class Cmpt>
+template<Foam::direction Row>
+inline void Foam::Tensor<Cmpt>::row(const Vector<Cmpt>& v)
+{
+    if (Row == 0)
+    {
+        this->v_[XX] = v.x(); this->v_[XY] = v.y(); this->v_[XZ] = v.z();
+    }
+    else if (Row == 1)
+    {
+        this->v_[YX] = v.x(); this->v_[YY] = v.y(); this->v_[YZ] = v.z();
+    }
+    else if (Row == 2)
+    {
+        this->v_[ZX] = v.x(); this->v_[ZY] = v.y(); this->v_[ZZ] = v.z();
+    }
+
+    static_assert(Row < 3, "Invalid row access");
+}
+
+
+template<class Cmpt>
+inline void Foam::Tensor<Cmpt>::col(const direction c, const Vector<Cmpt>& v)
+{
+    switch (c)
+    {
+        case 0: col<0>(v); break;
+        case 1: col<1>(v); break;
+        case 2: col<2>(v); break;
+        default:
+            FatalErrorInFunction
+                << "Invalid column access " << c << abort(FatalError);
+            break;
+    }
+}
+
+
+template<class Cmpt>
+inline void Foam::Tensor<Cmpt>::row(const direction r, const Vector<Cmpt>& v)
+{
+    switch (r)
+    {
+        case 0: row<0>(v); break;
+        case 1: row<1>(v); break;
+        case 2: row<2>(v); break;
+        default:
+            FatalErrorInFunction
+                << "Invalid row access " << r << abort(FatalError);
+            break;
+    }
+}
+
+
 template<class Cmpt>
 inline Foam::Vector<Cmpt> Foam::Tensor<Cmpt>::vectorComponent
 (
     const direction cmpt
 ) const
 {
-    switch (cmpt)
-    {
-        case 0:
-            return x();
-            break;
-        case 1:
-            return y();
-            break;
-        case 2:
-            return z();
-            break;
-    }
-
-    FatalErrorInFunction
-        << "Unhandled component " << cmpt << abort(FatalError);
-
-    return x();
+    return row(cmpt);
 }
 
 
+// * * * * * * * * * * * * * * * Member Operations * * * * * * * * * * * * * //
+
 template<class Cmpt>
 inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::T() const
 {
@@ -334,28 +481,73 @@ inline Foam::Tensor<Cmpt> Foam::Tensor<Cmpt>::T() const
 }
 
 
+template<class Cmpt>
+inline Foam::Tensor<Cmpt>
+Foam::Tensor<Cmpt>::inner(const Tensor<Cmpt>& t2) const
+{
+    const Tensor<Cmpt>& t1 = *this;
+
+    return Tensor<Cmpt>
+    (
+        t1.xx()*t2.xx() + t1.xy()*t2.yx() + t1.xz()*t2.zx(),
+        t1.xx()*t2.xy() + t1.xy()*t2.yy() + t1.xz()*t2.zy(),
+        t1.xx()*t2.xz() + t1.xy()*t2.yz() + t1.xz()*t2.zz(),
+
+        t1.yx()*t2.xx() + t1.yy()*t2.yx() + t1.yz()*t2.zx(),
+        t1.yx()*t2.xy() + t1.yy()*t2.yy() + t1.yz()*t2.zy(),
+        t1.yx()*t2.xz() + t1.yy()*t2.yz() + t1.yz()*t2.zz(),
+
+        t1.zx()*t2.xx() + t1.zy()*t2.yx() + t1.zz()*t2.zx(),
+        t1.zx()*t2.xy() + t1.zy()*t2.yy() + t1.zz()*t2.zy(),
+        t1.zx()*t2.xz() + t1.zy()*t2.yz() + t1.zz()*t2.zz()
+    );
+}
+
+
+template<class Cmpt>
+inline Foam::Tensor<Cmpt>
+Foam::Tensor<Cmpt>::innerT(const Tensor<Cmpt>& t2) const
+{
+    const Tensor<Cmpt>& t1 = *this;
+
+    return Tensor<Cmpt>
+    (
+        t1.xx()*t2.xx() + t1.xy()*t2.xy() + t1.xz()*t2.xz(),
+        t1.xx()*t2.yx() + t1.xy()*t2.yy() + t1.xz()*t2.yz(),
+        t1.xx()*t2.zx() + t1.xy()*t2.zy() + t1.xz()*t2.zz(),
+
+        t1.yx()*t2.xx() + t1.yy()*t2.xy() + t1.yz()*t2.xz(),
+        t1.yx()*t2.yx() + t1.yy()*t2.yy() + t1.yz()*t2.yz(),
+        t1.yx()*t2.zx() + t1.yy()*t2.zy() + t1.yz()*t2.zz(),
+
+        t1.zx()*t2.xx() + t1.zy()*t2.xy() + t1.zz()*t2.xz(),
+        t1.zx()*t2.yx() + t1.zy()*t2.yy() + t1.zz()*t2.yz(),
+        t1.zx()*t2.zx() + t1.zy()*t2.zy() + t1.zz()*t2.zz()
+    );
+}
+
+
+template<class Cmpt>
+inline Foam::Tensor<Cmpt>
+Foam::Tensor<Cmpt>::schur(const Tensor<Cmpt>& t2) const
+{
+    const Tensor<Cmpt>& t1 = *this;
+
+    return Tensor<Cmpt>
+    (
+        t1.xx()*t2.xx(), t1.xy()*t2.xy(), t1.xz()*t2.xz(),
+        t1.yx()*t2.yx(), t1.yy()*t2.yy(), t1.yz()*t2.yz(),
+        t1.zx()*t2.zx(), t1.zy()*t2.zy(), t1.zz()*t2.zz()
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Cmpt>
 inline void Foam::Tensor<Cmpt>::operator&=(const Tensor<Cmpt>& t)
 {
-    *this =
-    (
-        Tensor<Cmpt>
-        (
-            this->xx()*t.xx() + this->xy()*t.yx() + this->xz()*t.zx(),
-            this->xx()*t.xy() + this->xy()*t.yy() + this->xz()*t.zy(),
-            this->xx()*t.xz() + this->xy()*t.yz() + this->xz()*t.zz(),
-
-            this->yx()*t.xx() + this->yy()*t.yx() + this->yz()*t.zx(),
-            this->yx()*t.xy() + this->yy()*t.yy() + this->yz()*t.zy(),
-            this->yx()*t.xz() + this->yy()*t.yz() + this->yz()*t.zz(),
-
-            this->zx()*t.xx() + this->zy()*t.yx() + this->zz()*t.zx(),
-            this->zx()*t.xy() + this->zy()*t.yy() + this->zz()*t.zy(),
-            this->zx()*t.xz() + this->zy()*t.yz() + this->zz()*t.zz()
-        )
-    );
+    *this = this->inner(t);
 }
 
 
@@ -435,20 +627,7 @@ template<class Cmpt>
 inline typename innerProduct<Tensor<Cmpt>, Tensor<Cmpt>>::type
 operator&(const Tensor<Cmpt>& t1, const Tensor<Cmpt>& t2)
 {
-    return Tensor<Cmpt>
-    (
-        t1.xx()*t2.xx() + t1.xy()*t2.yx() + t1.xz()*t2.zx(),
-        t1.xx()*t2.xy() + t1.xy()*t2.yy() + t1.xz()*t2.zy(),
-        t1.xx()*t2.xz() + t1.xy()*t2.yz() + t1.xz()*t2.zz(),
-
-        t1.yx()*t2.xx() + t1.yy()*t2.yx() + t1.yz()*t2.zx(),
-        t1.yx()*t2.xy() + t1.yy()*t2.yy() + t1.yz()*t2.zy(),
-        t1.yx()*t2.xz() + t1.yy()*t2.yz() + t1.yz()*t2.zz(),
-
-        t1.zx()*t2.xx() + t1.zy()*t2.yx() + t1.zz()*t2.zx(),
-        t1.zx()*t2.xy() + t1.zy()*t2.yy() + t1.zz()*t2.zy(),
-        t1.zx()*t2.xz() + t1.zy()*t2.yz() + t1.zz()*t2.zz()
-    );
+    return t1.inner(t2);
 }
 
 
diff --git a/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H b/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
index 5ee3db5827..38de0da433 100644
--- a/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
+++ b/src/OpenFOAM/primitives/Tensor2D/Tensor2D.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -97,7 +97,7 @@ public:
         //- Construct given SphericalTensor2D
         inline Tensor2D(const SphericalTensor2D<Cmpt>&);
 
-        //- Construct given the two vectors
+        //- Construct given the two row vectors
         inline Tensor2D
         (
             const Vector2D<Cmpt>& x,
@@ -117,7 +117,7 @@ public:
 
     // Member Functions
 
-        // Access
+        // Component access
 
             inline const Cmpt& xx() const;
             inline const Cmpt& xy() const;
@@ -129,14 +129,41 @@ public:
             inline Cmpt& yx();
             inline Cmpt& yy();
 
-            // Access vector components.
 
+        // Column-vector access.
+
+            //- Extract vector for column 0
+            inline Vector2D<Cmpt> cx() const;
+
+            //- Extract vector for column 1
+            inline Vector2D<Cmpt> cy() const;
+
+
+        // Row-vector access.
+
+            //- Extract vector for row 0
             inline Vector2D<Cmpt> x() const;
+
+            //- Extract vector for row 1
             inline Vector2D<Cmpt> y() const;
 
+
+    // Tensor Operations
+
         //- Transpose
         inline Tensor2D<Cmpt> T() const;
 
+        //- Inner-product of this with another Tensor2D.
+        inline Tensor2D<Cmpt> inner(const Tensor2D<Cmpt>& t2) const;
+
+        //- Inner-product of this with transpose of another Tensor2D.
+        //  Primarily useful for coordinate transformations
+        //  (where transpose is the same as the inverse).
+        inline Tensor2D<Cmpt> innerT(const Tensor2D<Cmpt>& t2) const;
+
+        //- Schur-product of this with another Tensor2D.
+        inline Tensor2D<Cmpt> schur(const Tensor2D<Cmpt>& t2) const;
+
 
     // Member Operators
 
diff --git a/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H b/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
index 7f6be0d324..147155035f 100644
--- a/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
+++ b/src/OpenFOAM/primitives/Tensor2D/Tensor2DI.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -96,19 +96,6 @@ inline Foam::Tensor2D<Cmpt>::Tensor2D(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Cmpt>
-inline Foam::Vector2D<Cmpt> Foam::Tensor2D<Cmpt>::x() const
-{
-    return Vector2D<Cmpt>(this->v_[XX], this->v_[XY]);
-}
-
-template<class Cmpt>
-inline Foam::Vector2D<Cmpt> Foam::Tensor2D<Cmpt>::y() const
-{
-    return Vector2D<Cmpt>(this->v_[YX], this->v_[YY]);
-}
-
-
 template<class Cmpt>
 inline const Cmpt& Foam::Tensor2D<Cmpt>::xx() const
 {
@@ -159,6 +146,34 @@ inline Cmpt& Foam::Tensor2D<Cmpt>::yy()
 }
 
 
+template<class Cmpt>
+inline Foam::Vector2D<Cmpt> Foam::Tensor2D<Cmpt>::x() const
+{
+    return Vector2D<Cmpt>(this->v_[XX], this->v_[XY]);
+}
+
+template<class Cmpt>
+inline Foam::Vector2D<Cmpt> Foam::Tensor2D<Cmpt>::y() const
+{
+    return Vector2D<Cmpt>(this->v_[YX], this->v_[YY]);
+}
+
+
+template<class Cmpt>
+inline Foam::Vector2D<Cmpt> Foam::Tensor2D<Cmpt>::cx() const
+{
+    return Vector2D<Cmpt>(this->v_[XX], this->v_[YX]);
+}
+
+template<class Cmpt>
+inline Foam::Vector2D<Cmpt> Foam::Tensor2D<Cmpt>::cy() const
+{
+    return Vector2D<Cmpt>(this->v_[XY], this->v_[YY]);
+}
+
+
+// * * * * * * * * * * * * * * * Member Operations * * * * * * * * * * * * * //
+
 template<class Cmpt>
 inline Foam::Tensor2D<Cmpt> Foam::Tensor2D<Cmpt>::T() const
 {
@@ -170,6 +185,54 @@ inline Foam::Tensor2D<Cmpt> Foam::Tensor2D<Cmpt>::T() const
 }
 
 
+template<class Cmpt>
+inline Foam::Tensor2D<Cmpt>
+Foam::Tensor2D<Cmpt>::inner(const Tensor2D<Cmpt>& t2) const
+{
+    const Tensor2D<Cmpt>& t1 = *this;
+
+    return Tensor2D<Cmpt>
+    (
+        t1.xx()*t2.xx() + t1.xy()*t2.yx(),
+        t1.xx()*t2.xy() + t1.xy()*t2.yy(),
+
+        t1.yx()*t2.xx() + t1.yy()*t2.yx(),
+        t1.yx()*t2.xy() + t1.yy()*t2.yy()
+    );
+}
+
+
+template<class Cmpt>
+inline Foam::Tensor2D<Cmpt>
+Foam::Tensor2D<Cmpt>::innerT(const Tensor2D<Cmpt>& t2) const
+{
+    const Tensor2D<Cmpt>& t1 = *this;
+
+    return Tensor2D<Cmpt>
+    (
+        t1.xx()*t2.xx() + t1.xy()*t2.xy(),
+        t1.xx()*t2.yx() + t1.xy()*t2.yy(),
+
+        t1.yx()*t2.xx() + t1.yy()*t2.xy(),
+        t1.yx()*t2.yx() + t1.yy()*t2.yy()
+    );
+}
+
+
+template<class Cmpt>
+inline Foam::Tensor2D<Cmpt>
+Foam::Tensor2D<Cmpt>::schur(const Tensor2D<Cmpt>& t2) const
+{
+    const Tensor2D<Cmpt>& t1 = *this;
+
+    return Tensor2D<Cmpt>
+    (
+        t1.xx()*t2.xx(), t1.xy()*t2.xy(),
+        t1.yx()*t2.yx(), t1.yy()*t2.yy()
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class Cmpt>
@@ -200,16 +263,10 @@ template<class Cmpt>
 inline typename innerProduct<Tensor2D<Cmpt>, Tensor2D<Cmpt>>::type
 operator&(const Tensor2D<Cmpt>& t1, const Tensor2D<Cmpt>& t2)
 {
-    return Tensor2D<Cmpt>
-    (
-        t1.xx()*t2.xx() + t1.xy()*t2.yx(),
-        t1.xx()*t2.xy() + t1.xy()*t2.yy(),
-
-        t1.yx()*t2.xx() + t1.yy()*t2.yx(),
-        t1.yx()*t2.xy() + t1.yy()*t2.yy()
-    );
+    return t1.inner(t2);
 }
 
+
 //- Inner-product between a tensor and a vector
 template<class Cmpt>
 inline typename innerProduct<Tensor2D<Cmpt>, Vector2D<Cmpt>>::type
diff --git a/src/OpenFOAM/primitives/transform/transform.H b/src/OpenFOAM/primitives/transform/transform.H
index 9fd1318acc..cb72b2ba18 100644
--- a/src/OpenFOAM/primitives/transform/transform.H
+++ b/src/OpenFOAM/primitives/transform/transform.H
@@ -142,31 +142,68 @@ inline tensor Ra(const vector& a, const scalar omega)
 }
 
 
-inline label transform(const tensor&, const bool i)
+//- No-op rotational transform of a bool
+inline bool transform(const tensor&, const bool b)
 {
-    return i;
+    return b;
 }
 
 
+//- No-op inverse rotational transform of a bool
+inline bool invTransform(const tensor&, const bool b)
+{
+    return b;
+}
+
+
+//- No-op rotational transform of a label
 inline label transform(const tensor&, const label i)
 {
     return i;
 }
 
 
+//- No-op inverse rotational transform of a label
+inline label invTransform(const tensor&, const label i)
+{
+    return i;
+}
+
+
+//- No-op rotational transform of a label
 inline scalar transform(const tensor&, const scalar s)
 {
     return s;
 }
 
 
-template<class Cmpt>
-inline Vector<Cmpt> transform(const tensor& tt, const Vector<Cmpt>& v)
+//- No-op inverse rotational transform of a scalar
+inline scalar invTransform(const tensor&, const scalar s)
 {
-    return tt & v;
+    return s;
 }
 
 
+//- Use rotational tensor to transform a vector.
+//  Same as (rot & v)
+template<class Cmpt>
+inline Vector<Cmpt> transform(const tensor& tt, const Vector<Cmpt>& v)
+{
+    return (tt & v);
+}
+
+
+//- Use rotational tensor to inverse transform a vector.
+//  Same as (v & rot)
+template<class Cmpt>
+inline Vector<Cmpt> invTransform(const tensor& tt, const Vector<Cmpt>& v)
+{
+    return (v & tt);
+}
+
+
+//- Use rotational tensor to transform a tensor.
+//  Same as (rot & input & rot.T())
 template<class Cmpt>
 inline Tensor<Cmpt> transform(const tensor& tt, const Tensor<Cmpt>& t)
 {
@@ -211,6 +248,53 @@ inline Tensor<Cmpt> transform(const tensor& tt, const Tensor<Cmpt>& t)
 }
 
 
+//- Use rotational tensor to inverse transform a tensor.
+//  Same as (rot.T() & input & rot)
+template<class Cmpt>
+inline Tensor<Cmpt> invTransform(const tensor& tt, const Tensor<Cmpt>& t)
+{
+    return Tensor<Cmpt>
+    (
+        (tt.xx()*t.xx() + tt.yx()*t.yx() + tt.zx()*t.zx())*tt.xx()
+      + (tt.xx()*t.xy() + tt.yx()*t.yy() + tt.zx()*t.zy())*tt.yx()
+      + (tt.xx()*t.xz() + tt.yx()*t.yz() + tt.zx()*t.zz())*tt.zx(),
+
+        (tt.xx()*t.xx() + tt.yx()*t.yx() + tt.zx()*t.zx())*tt.xy()
+      + (tt.xx()*t.xy() + tt.yx()*t.yy() + tt.zx()*t.zy())*tt.yy()
+      + (tt.xx()*t.xz() + tt.yx()*t.yz() + tt.zx()*t.zz())*tt.zy(),
+
+        (tt.xx()*t.xx() + tt.yx()*t.yx() + tt.zx()*t.zx())*tt.xz()
+      + (tt.xx()*t.xy() + tt.yx()*t.yy() + tt.zx()*t.zy())*tt.yz()
+      + (tt.xx()*t.xz() + tt.yx()*t.yz() + tt.zx()*t.zz())*tt.zz(),
+
+        (tt.xy()*t.xx() + tt.yy()*t.yx() + tt.zy()*t.zx())*tt.xx()
+      + (tt.xy()*t.xy() + tt.yy()*t.yy() + tt.zy()*t.zy())*tt.yx()
+      + (tt.xy()*t.xz() + tt.yy()*t.yz() + tt.zy()*t.zz())*tt.zx(),
+
+        (tt.xy()*t.xx() + tt.yy()*t.yx() + tt.zy()*t.zx())*tt.xy()
+      + (tt.xy()*t.xy() + tt.yy()*t.yy() + tt.zy()*t.zy())*tt.yy()
+      + (tt.xy()*t.xz() + tt.yy()*t.yz() + tt.zy()*t.zz())*tt.zy(),
+
+        (tt.xy()*t.xx() + tt.yy()*t.yx() + tt.zy()*t.zx())*tt.xz()
+      + (tt.xy()*t.xy() + tt.yy()*t.yy() + tt.zy()*t.zy())*tt.yz()
+      + (tt.xy()*t.xz() + tt.yy()*t.yz() + tt.zy()*t.zz())*tt.zz(),
+
+        (tt.xz()*t.xx() + tt.yz()*t.yx() + tt.zz()*t.zx())*tt.xx()
+      + (tt.xz()*t.xy() + tt.yz()*t.yy() + tt.zz()*t.zy())*tt.yx()
+      + (tt.xz()*t.xz() + tt.yz()*t.yz() + tt.zz()*t.zz())*tt.zx(),
+
+        (tt.xz()*t.xx() + tt.yz()*t.yx() + tt.zz()*t.zx())*tt.xy()
+      + (tt.xz()*t.xy() + tt.yz()*t.yy() + tt.zz()*t.zy())*tt.yy()
+      + (tt.xz()*t.xz() + tt.yz()*t.yz() + tt.zz()*t.zz())*tt.zy(),
+
+        (tt.xz()*t.xx() + tt.yz()*t.yx() + tt.zz()*t.zx())*tt.xz()
+      + (tt.xz()*t.xy() + tt.yz()*t.yy() + tt.zz()*t.zy())*tt.yz()
+      + (tt.xz()*t.xz() + tt.yz()*t.yz() + tt.zz()*t.zz())*tt.zz()
+    );
+}
+
+
+//- Use rotational tensor to transform a spherical tensor (no-op).
 template<class Cmpt>
 inline SphericalTensor<Cmpt> transform
 (
@@ -222,6 +306,20 @@ inline SphericalTensor<Cmpt> transform
 }
 
 
+//- Use rotational tensor to inverse transform a spherical tensor (no-op).
+template<class Cmpt>
+inline SphericalTensor<Cmpt> invTransform
+(
+    const tensor& tt,
+    const SphericalTensor<Cmpt>& st
+)
+{
+    return st;
+}
+
+
+//- Use rotational tensor to transform a symmetrical tensor.
+//  Same as (rot & input & rot.T())
 template<class Cmpt>
 inline SymmTensor<Cmpt> transform(const tensor& tt, const SymmTensor<Cmpt>& st)
 {
@@ -254,6 +352,41 @@ inline SymmTensor<Cmpt> transform(const tensor& tt, const SymmTensor<Cmpt>& st)
 }
 
 
+//- Use rotational tensor to inverse transform a symmetrical tensor.
+//  Same as (rot.T() & input & rot)
+template<class Cmpt>
+inline SymmTensor<Cmpt>
+invTransform(const tensor& tt, const SymmTensor<Cmpt>& st)
+{
+    return SymmTensor<Cmpt>
+    (
+        (tt.xx()*st.xx() + tt.yx()*st.xy() + tt.zx()*st.xz())*tt.xx()
+      + (tt.xx()*st.xy() + tt.yx()*st.yy() + tt.zx()*st.yz())*tt.yx()
+      + (tt.xx()*st.xz() + tt.yx()*st.yz() + tt.zx()*st.zz())*tt.zx(),
+
+        (tt.xx()*st.xx() + tt.yx()*st.xy() + tt.zx()*st.xz())*tt.xy()
+      + (tt.xx()*st.xy() + tt.yx()*st.yy() + tt.zx()*st.yz())*tt.yy()
+      + (tt.xx()*st.xz() + tt.yx()*st.yz() + tt.zx()*st.zz())*tt.zy(),
+
+        (tt.xx()*st.xx() + tt.yx()*st.xy() + tt.zx()*st.xz())*tt.xz()
+      + (tt.xx()*st.xy() + tt.yx()*st.yy() + tt.zx()*st.yz())*tt.yz()
+      + (tt.xx()*st.xz() + tt.yx()*st.yz() + tt.zx()*st.zz())*tt.zz(),
+
+        (tt.xy()*st.xx() + tt.yy()*st.xy() + tt.zy()*st.xz())*tt.xy()
+      + (tt.xy()*st.xy() + tt.yy()*st.yy() + tt.zy()*st.yz())*tt.yy()
+      + (tt.xy()*st.xz() + tt.yy()*st.yz() + tt.zy()*st.zz())*tt.zy(),
+
+        (tt.xy()*st.xx() + tt.yy()*st.xy() + tt.zy()*st.xz())*tt.xz()
+      + (tt.xy()*st.xy() + tt.yy()*st.yy() + tt.zy()*st.yz())*tt.yz()
+      + (tt.xy()*st.xz() + tt.yy()*st.yz() + tt.zy()*st.zz())*tt.zz(),
+
+        (tt.xz()*st.xx() + tt.yz()*st.xy() + tt.zz()*st.xz())*tt.xz()
+      + (tt.xz()*st.xy() + tt.yz()*st.yy() + tt.zz()*st.yz())*tt.yz()
+      + (tt.xz()*st.xz() + tt.yz()*st.yz() + tt.zz()*st.zz())*tt.zz()
+    );
+}
+
+
 template<class Type1, class Type2>
 inline Type1 transformMask(const Type2& t)
 {
@@ -292,16 +425,16 @@ inline scalar pseudoAngle
     const vector& vec
 )
 {
-    scalar cos = vec & e0;
-    scalar sin = vec & e1;
+    const scalar cos_angle = vec & e0;
+    const scalar sin_angle = vec & e1;
 
-    if (sin < -SMALL)
+    if (sin_angle < -SMALL)
     {
-        return (3.0 + cos)*constant::mathematical::piByTwo;
+        return (3.0 + cos_angle)*constant::mathematical::piByTwo;
     }
     else
     {
-        return (1.0 - cos)*constant::mathematical::piByTwo;
+        return (1.0 - cos_angle)*constant::mathematical::piByTwo;
     }
 }
 
diff --git a/src/OpenFOAM/primitives/triad/triad.C b/src/OpenFOAM/primitives/triad/triad.C
index 0a54dc4aca..ec800e9f3d 100644
--- a/src/OpenFOAM/primitives/triad/triad.C
+++ b/src/OpenFOAM/primitives/triad/triad.C
@@ -95,14 +95,6 @@ Foam::triad::triad(const quaternion& q)
 }
 
 
-Foam::triad::triad(const tensor& t)
-{
-    x() = t.x();
-    y() = t.y();
-    z() = t.z();
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::triad::orthogonalize()
@@ -379,16 +371,6 @@ Foam::triad::operator Foam::quaternion() const
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-void Foam::triad::operator=(const tensor& t)
-{
-    x() = t.x();
-    y() = t.y();
-    z() = t.z();
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 Foam::scalar Foam::diff(const triad& A, const triad& B)
diff --git a/src/OpenFOAM/primitives/triad/triad.H b/src/OpenFOAM/primitives/triad/triad.H
index 30eaba81c1..9a1acd1ede 100644
--- a/src/OpenFOAM/primitives/triad/triad.H
+++ b/src/OpenFOAM/primitives/triad/triad.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,7 +25,8 @@ Class
     Foam::triad
 
 Description
-    Representation of a 3D Cartesian coordinate system as a Vector of vectors.
+    Representation of a 3D Cartesian coordinate system as a Vector of
+    row vectors.
 
 See also
     Foam::quaternion
@@ -48,12 +49,10 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward declarations
+class triad;
 class Istream;
 class Ostream;
-
-// Forward declaration of friend functions and operators
-class triad;
 Istream& operator>>(Istream&, triad&);
 Ostream& operator<<(Ostream&, const triad&);
 
@@ -78,18 +77,18 @@ public:
         //- Construct from components
         inline triad(const Vector<vector>& vv);
 
-        //- Construct from coordinate axes
+        //- Construct from coordinate axes (row vectors)
         inline triad(const vector& x, const vector& y, const vector& z);
 
+        //- Construct from a tensor
+        inline triad(const tensor& t);
+
         //- Construct from a primary axis with the other two unset
         inline triad(const vector& pa);
 
         //- Construct from a quaternion
         triad(const quaternion& q);
 
-        //- Construct from a tensor
-        triad(const tensor& t);
-
         //- Construct from Istream
         inline triad(Istream&);
 
@@ -130,12 +129,27 @@ public:
         //- Convert to a quaternion
         operator quaternion() const;
 
+        //- Return transpose
+        inline triad T() const;
+
+
+    // Column-vector access.
+
+        //- Extract vector for column 0
+        inline vector cx() const;
+
+        //- Extract vector for column 1
+        inline vector cy() const;
+
+        //- Extract vector for column 2
+        inline vector cz() const;
+
 
     // Member Operators
 
-        inline void operator=(const Vector<vector>&);
+        inline void operator=(const Vector<vector>& vv);
 
-        void operator=(const tensor& t);
+        inline void operator=(const tensor& t);
 
         //- Add the triad t2 to this triad
         //  without normalizing or orthogonalizing
diff --git a/src/OpenFOAM/primitives/triad/triadI.H b/src/OpenFOAM/primitives/triad/triadI.H
index d259258178..249f5f9157 100644
--- a/src/OpenFOAM/primitives/triad/triadI.H
+++ b/src/OpenFOAM/primitives/triad/triadI.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,6 +43,12 @@ inline Foam::triad::triad(const vector& x, const vector& y, const vector& z)
 {}
 
 
+inline Foam::triad::triad(const tensor& t)
+:
+    Vector<vector>(t.x(), t.y(), t.z())
+{}
+
+
 inline Foam::triad::triad(const vector& pa)
 {
     operator=(triad::unset);
@@ -116,6 +122,30 @@ inline void Foam::triad::normalize()
 }
 
 
+inline Foam::vector Foam::triad::cx() const
+{
+    return vector(x().x(), y().x(), z().x());
+}
+
+
+inline Foam::vector Foam::triad::cy() const
+{
+    return vector(x().y(), y().y(), z().y());
+}
+
+
+inline Foam::vector Foam::triad::cz() const
+{
+    return vector(x().z(), y().z(), z().z());
+}
+
+
+inline Foam::triad Foam::triad::T() const
+{
+    return triad(cx(), cy(), cz());
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 inline void Foam::triad::operator=(const Vector<vector>& vv)
@@ -124,6 +154,14 @@ inline void Foam::triad::operator=(const Vector<vector>& vv)
 }
 
 
+inline void Foam::triad::operator=(const tensor& t)
+{
+    x() = t.x();
+    y() = t.y();
+    z() = t.z();
+}
+
+
 // * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
 
 inline Foam::Istream& Foam::operator>>(Istream& is, triad& t)
diff --git a/src/conversion/ensight/output/ensightOutputTemplates.C b/src/conversion/ensight/output/ensightOutputTemplates.C
index a29fa3ea58..3ffb066149 100644
--- a/src/conversion/ensight/output/ensightOutputTemplates.C
+++ b/src/conversion/ensight/output/ensightOutputTemplates.C
@@ -261,7 +261,7 @@ bool Foam::ensightOutput::writeField
         // flat boundary field
         // as per volPointInterpolation::flatBoundaryField()
 
-        Field<Type> flat(mesh.nFaces() - mesh.nInternalFaces());
+        Field<Type> flat(mesh.nBoundaryFaces());
 
         const fvBoundaryMesh& bm = mesh.boundary();
         forAll(vf.boundaryField(), patchI)
diff --git a/src/conversion/polyDualMesh/polyDualMesh.C b/src/conversion/polyDualMesh/polyDualMesh.C
index 9b325c9415..382ae4d445 100644
--- a/src/conversion/polyDualMesh/polyDualMesh.C
+++ b/src/conversion/polyDualMesh/polyDualMesh.C
@@ -1407,7 +1407,7 @@ Foam::polyDualMesh::polyDualMesh
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        labelList(mesh.nFaces() - mesh.nInternalFaces())
+        labelList(mesh.nBoundaryFaces(), -1)
     )
 {
     calcDual(mesh, featureEdges, featurePoints);
@@ -1446,7 +1446,7 @@ Foam::polyDualMesh::polyDualMesh
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        labelList(mesh.nFaces() - mesh.nInternalFaces(), -1)
+        labelList(mesh.nBoundaryFaces(), -1)
     )
 {
     labelList featureEdges, featurePoints;
@@ -1470,7 +1470,7 @@ void Foam::polyDualMesh::calcFeatures
         SubList<face>
         (
             mesh.faces(),
-            mesh.nFaces() - mesh.nInternalFaces(),
+            mesh.nBoundaryFaces(),
             mesh.nInternalFaces()
         ),
         mesh.points()
diff --git a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
index 25dc0c9687..266d2c3dbf 100644
--- a/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
+++ b/src/dynamicFvMesh/dynamicRefineFvMesh/dynamicRefineFvMesh.C
@@ -86,7 +86,7 @@ void Foam::dynamicRefineFvMesh::calculateProtectedCells
     unrefineableCell = protectedCell_;
 
     // Get neighbouring cell level
-    labelList neiLevel(nFaces()-nInternalFaces());
+    labelList neiLevel(nBoundaryFaces());
 
     for (label facei = nInternalFaces(); facei < nFaces(); facei++)
     {
diff --git a/src/dynamicMesh/boundaryMesh/boundaryMesh.C b/src/dynamicMesh/boundaryMesh/boundaryMesh.C
index a26a431330..1d3a1ec380 100644
--- a/src/dynamicMesh/boundaryMesh/boundaryMesh.C
+++ b/src/dynamicMesh/boundaryMesh/boundaryMesh.C
@@ -476,7 +476,7 @@ void Foam::boundaryMesh::read(const polyMesh& mesh)
     patches_.setSize(mesh.boundaryMesh().size());
 
     // Number of boundary faces
-    label nBFaces = mesh.nFaces() - mesh.nInternalFaces();
+    const label nBFaces = mesh.nBoundaryFaces();
 
     faceList bFaces(nBFaces);
 
@@ -976,7 +976,7 @@ Foam::labelList Foam::boundaryMesh::getNearest
     // Search nearest triangle centre for every polyMesh boundary face
     //
 
-    labelList nearestBFacei(pMesh.nFaces() - pMesh.nInternalFaces());
+    labelList nearestBFacei(pMesh.nBoundaryFaces());
 
     treeBoundBox tightest;
 
@@ -1168,7 +1168,7 @@ void Foam::boundaryMesh::patchify
     label meshFacei = newMesh.nInternalFaces();
 
     // First patch gets all non-coupled faces
-    label facesToBeDone = newMesh.nFaces() - newMesh.nInternalFaces();
+    label facesToBeDone = newMesh.nBoundaryFaces();
 
     forAll(patches_, bPatchi)
     {
@@ -1262,9 +1262,7 @@ void Foam::boundaryMesh::patchify
         List<DynamicList<label>> patchFaces(nNewPatches);
 
         // Give reasonable estimate for size of patches
-        label nAvgFaces =
-            (newMesh.nFaces() - newMesh.nInternalFaces())
-          / nNewPatches;
+        label nAvgFaces = newMesh.nBoundaryFaces() / nNewPatches;
 
         forAll(patchFaces, newPatchi)
         {
diff --git a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
index d6c18d1cb8..02ef7638ea 100644
--- a/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
+++ b/src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C
@@ -500,7 +500,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::deleteProcPatches
 
     // New patchID per boundary faces to be repatched. Is -1 (no change)
     // or new patchID
-    labelList newPatchID(mesh_.nFaces() - mesh_.nInternalFaces(), -1);
+    labelList newPatchID(mesh_.nBoundaryFaces(), -1);
 
     label nProcPatches = 0;
 
@@ -791,7 +791,7 @@ void Foam::fvMeshDistribute::getCouplingData
     // Construct the coupling information for all (boundary) faces and
     // points
 
-    label nBnd = mesh_.nFaces() - mesh_.nInternalFaces();
+    const label nBnd = mesh_.nBoundaryFaces();
     sourceFace.setSize(nBnd);
     sourceProc.setSize(nBnd);
     sourcePatch.setSize(nBnd);
@@ -998,10 +998,10 @@ void Foam::fvMeshDistribute::subsetCouplingData
     labelList& subPointMaster
 )
 {
-    subFace.setSize(mesh.nFaces() - mesh.nInternalFaces());
-    subProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
-    subPatch.setSize(mesh.nFaces() - mesh.nInternalFaces());
-    subNewNbrProc.setSize(mesh.nFaces() - mesh.nInternalFaces());
+    subFace.setSize(mesh.nBoundaryFaces());
+    subProc.setSize(mesh.nBoundaryFaces());
+    subPatch.setSize(mesh.nBoundaryFaces());
+    subNewNbrProc.setSize(mesh.nBoundaryFaces());
 
     forAll(subFace, newBFacei)
     {
@@ -1132,7 +1132,7 @@ Foam::labelList Foam::fvMeshDistribute::mapBoundaryData
     const labelList& boundaryData1  // on added mesh
 )
 {
-    labelList newBoundaryData(mesh.nFaces() - mesh.nInternalFaces());
+    labelList newBoundaryData(mesh.nBoundaryFaces());
 
     forAll(boundaryData0, oldBFacei)
     {
@@ -1932,7 +1932,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
             SubList<label>
             (
                 repatchMap().reverseFaceMap(),
-                mesh_.nFaces() - mesh_.nInternalFaces(),
+                mesh_.nBoundaryFaces(),
                 mesh_.nInternalFaces()
             )
           - mesh_.nInternalFaces()
@@ -2685,7 +2685,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
             {
                 // Find out if any faces of domain mesh were flipped (boundary
                 // faces becoming internal)
-                label nBnd = domainMesh.nFaces()-domainMesh.nInternalFaces();
+                const label nBnd = domainMesh.nBoundaryFaces();
                 flippedAddedFaces.resize(nBnd/4);
 
                 for
diff --git a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
index edb8209d6f..97b23b8bf7 100644
--- a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
+++ b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
@@ -139,7 +139,7 @@ void Foam::cellCuts::syncProc()
     syncTools::syncEdgeList(mesh(), edgeWeight_, maxEqOp<scalar>(), -GREAT);
 
     {
-        const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+        const label nBnd = mesh().nBoundaryFaces();
 
         // Convert faceSplitCut into face-local data: vertex and edge w.r.t.
         // vertex 0: (since this is same on both sides)
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 2a0b294eba..46500f1399 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
@@ -385,7 +385,7 @@ bool Foam::polyMeshGeometry::checkFaceDotProduct
     const scalar severeNonorthogonalityThreshold = ::cos(degToRad(orthWarn));
 
     // Calculate coupled cell centre
-    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
+    pointField neiCc(mesh.nBoundaryFaces());
 
     for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
     {
@@ -745,7 +745,7 @@ bool Foam::polyMeshGeometry::checkFaceTets
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     // Calculate coupled cell centre
-    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
+    pointField neiCc(mesh.nBoundaryFaces());
 
     for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
     {
@@ -1187,7 +1187,7 @@ bool Foam::polyMeshGeometry::checkFaceWeights
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     // Calculate coupled cell centre
-    pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
+    pointField neiCc(mesh.nBoundaryFaces());
 
     for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
     {
@@ -1340,7 +1340,7 @@ bool Foam::polyMeshGeometry::checkVolRatio
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     // Calculate coupled cell vol
-    scalarField neiVols(mesh.nFaces()-mesh.nInternalFaces());
+    scalarField neiVols(mesh.nBoundaryFaces());
 
     for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
     {
@@ -1671,7 +1671,7 @@ bool Foam::polyMeshGeometry::checkFaceTwist
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
     // Calculate coupled cell centre
-    pointField neiCc(mesh.nFaces()-mesh.nInternalFaces());
+    pointField neiCc(mesh.nBoundaryFaces());
 
     for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
     {
diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
index be4b283afd..a7135b9f4f 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
+++ b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
@@ -908,7 +908,7 @@ void Foam::faceCoupleInfo::findPerfectMatchingFaces
             mesh0.faces(),
             mesh0.points(),
             mesh0.nInternalFaces(),
-            mesh0.nFaces() - mesh0.nInternalFaces()
+            mesh0.nBoundaryFaces()
         )
     );
 
@@ -919,7 +919,7 @@ void Foam::faceCoupleInfo::findPerfectMatchingFaces
             mesh1.faces(),
             mesh1.points(),
             mesh1.nInternalFaces(),
-            mesh1.nFaces() - mesh1.nInternalFaces()
+            mesh1.nBoundaryFaces()
         )
     );
 
@@ -984,11 +984,10 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
 )
 {
     // Construct octree from all mesh0 boundary faces
-    labelList bndFaces(mesh0.nFaces()-mesh0.nInternalFaces());
-    forAll(bndFaces, i)
-    {
-        bndFaces[i] = mesh0.nInternalFaces() + i;
-    }
+    labelList bndFaces
+    (
+        identity(mesh0.nBoundaryFaces(), mesh0.nInternalFaces())
+    );
 
     treeBoundBox overallBb(mesh0.points());
 
@@ -1016,8 +1015,8 @@ void Foam::faceCoupleInfo::findSlavesCoveringMaster
 
     // Search nearest face for every mesh1 boundary face.
 
-    labelHashSet mesh0Set(mesh0.nFaces() - mesh0.nInternalFaces());
-    labelHashSet mesh1Set(mesh1.nFaces() - mesh1.nInternalFaces());
+    labelHashSet mesh0Set(mesh0.nBoundaryFaces());
+    labelHashSet mesh1Set(mesh1.nBoundaryFaces());
 
     for
     (
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
index 44c3ca0afc..34ad10fde8 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C
@@ -384,7 +384,7 @@ Foam::labelListList Foam::combineFaces::getMergeSets
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
     // Pick up all cells on boundary
-    labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelHashSet boundaryCells(mesh_.nBoundaryFaces());
 
     forAll(patches, patchi)
     {
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
index 0608d45f7b..dae24bc067 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
@@ -462,7 +462,7 @@ Foam::scalarField Foam::edgeCollapser::calcTargetFaceSizes() const
     const labelList& cellOwner = mesh_.faceOwner();
     const labelList& cellNeighbour = mesh_.faceNeighbour();
 
-    const label nBoundaryFaces = mesh_.nFaces() - mesh_.nInternalFaces();
+    const label nBoundaryFaces = mesh_.nBoundaryFaces();
 
     // Calculate face size from cell volumes for internal faces
     for (label intFacei = 0; intFacei < mesh_.nInternalFaces(); ++intFacei)
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
index 8d1b96c318..fa3797e2f7 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/hexRef8/hexRef8.C
@@ -1608,7 +1608,7 @@ Foam::label Foam::hexRef8::faceConsistentRefinement
 
     // Coupled faces. Swap owner level to get neighbouring cell level.
     // (only boundary faces of neiLevel used)
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
 
     forAll(neiLevel, i)
     {
@@ -1685,7 +1685,7 @@ void Foam::hexRef8::checkWantedRefinementLevels
 
     // Coupled faces. Swap owner level to get neighbouring cell level.
     // (only boundary faces of neiLevel used)
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
 
     forAll(neiLevel, i)
     {
@@ -2696,9 +2696,9 @@ Foam::labelList Foam::hexRef8::consistentSlowRefinement
         // Coupled faces. Swap owner level to get neighbouring cell level.
         // (only boundary faces of neiLevel used)
 
-        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-        labelList neiCount(mesh_.nFaces()-mesh_.nInternalFaces());
-        labelList neiRefCount(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList neiLevel(mesh_.nBoundaryFaces());
+        labelList neiCount(mesh_.nBoundaryFaces());
+        labelList neiRefCount(mesh_.nBoundaryFaces());
 
         forAll(neiLevel, i)
         {
@@ -3471,7 +3471,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
     //  refinining and subsetting)
 
     {
-        labelList newNeiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList newNeiLevel(mesh_.nBoundaryFaces());
 
         forAll(newNeiLevel, i)
         {
@@ -3528,7 +3528,7 @@ Foam::labelListList Foam::hexRef8::setRefinement
         // above
         pointField bFaceMids
         (
-            mesh_.nFaces()-mesh_.nInternalFaces(),
+            mesh_.nBoundaryFaces(),
             point(-GREAT, -GREAT, -GREAT)
         );
 
@@ -4562,7 +4562,7 @@ void Foam::hexRef8::checkMesh() const
     // otherwise mesh redistribution might cause multiple faces between two
     // cells
     {
-        labelList nei(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList nei(mesh_.nBoundaryFaces());
         forAll(nei, i)
         {
             nei[i] = mesh_.faceOwner()[i+mesh_.nInternalFaces()];
@@ -4616,7 +4616,7 @@ void Foam::hexRef8::checkMesh() const
     // ~~~~~~~~~~~~~~~~~
 
     {
-        scalarField neiFaceAreas(mesh_.nFaces()-mesh_.nInternalFaces());
+        scalarField neiFaceAreas(mesh_.nBoundaryFaces());
         forAll(neiFaceAreas, i)
         {
             neiFaceAreas[i] = mag(mesh_.faceAreas()[i+mesh_.nInternalFaces()]);
@@ -4658,7 +4658,7 @@ void Foam::hexRef8::checkMesh() const
     // Check number of points on faces.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     {
-        labelList nVerts(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList nVerts(mesh_.nBoundaryFaces());
 
         forAll(nVerts, i)
         {
@@ -4700,7 +4700,7 @@ void Foam::hexRef8::checkMesh() const
     // ~~~~~~~~~~~~~~~~~~~~
     {
         // Anchor points.
-        pointField anchorPoints(mesh_.nFaces()-mesh_.nInternalFaces());
+        pointField anchorPoints(mesh_.nBoundaryFaces());
 
         forAll(anchorPoints, i)
         {
@@ -4808,7 +4808,7 @@ void Foam::hexRef8::checkRefinementLevels
         }
 
         // Coupled faces. Get neighbouring value
-        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList neiLevel(mesh_.nBoundaryFaces());
 
         forAll(neiLevel, i)
         {
@@ -5379,7 +5379,7 @@ Foam::labelList Foam::hexRef8::consistentUnrefinement
 
 
         // Coupled faces. Swap owner level to get neighbouring cell level.
-        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList neiLevel(mesh_.nBoundaryFaces());
 
         forAll(neiLevel, i)
         {
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C
index 38269d63a6..2e6244fe37 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C
@@ -133,7 +133,7 @@ Foam::labelList Foam::removeCells::getExposedFaces
         SubList<label> bndValues
         (
             nCellsUsingFace,
-            mesh_.nFaces()-mesh_.nInternalFaces(),
+            mesh_.nBoundaryFaces(),
             mesh_.nInternalFaces()
         );
 
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
index 4cd0f0e0ae..30c1382e82 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/removePoints.C
@@ -638,7 +638,7 @@ void Foam::removePoints::getUnrefimentSet
         // restoring. Note that this is over all saved faces, not just over
         // the ones in undoFaces.
 
-        boolListList faceVertexRestore(mesh_.nFaces()-mesh_.nInternalFaces());
+        boolListList faceVertexRestore(mesh_.nBoundaryFaces());
 
         // Populate with my local points-to-restore.
         forAll(savedFaces_, saveI)
diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
index dea8299a49..dde8008349 100644
--- a/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
+++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/tetDecomposer.C
@@ -222,7 +222,7 @@ void Foam::tetDecomposer::setRefinement
             }
         }
 
-        boolList neiDecomposeCell(mesh_.nFaces()-mesh_.nInternalFaces());
+        boolList neiDecomposeCell(mesh_.nBoundaryFaces());
         forAll(neiDecomposeCell, bFacei)
         {
             label facei = mesh_.nInternalFaces()+bFacei;
diff --git a/src/finiteArea/faMesh/faMeshMapper/faAreaMapper.C b/src/finiteArea/faMesh/faMeshMapper/faAreaMapper.C
index de29f2e4a0..0b772daad9 100644
--- a/src/finiteArea/faMesh/faMeshMapper/faAreaMapper.C
+++ b/src/finiteArea/faMesh/faMeshMapper/faAreaMapper.C
@@ -63,18 +63,10 @@ void Foam::faAreaMapper::calcAddressing() const
 
     // Prepare a list of new face labels and (preliminary) addressing
     // Note: dimensioned to number of boundary faces of polyMesh
-    newFaceLabelsPtr_ = new labelList
-    (
-        mesh_().nFaces() - mesh_().nInternalFaces(),
-        -1
-    );
+    newFaceLabelsPtr_ = new labelList(mesh_().nBoundaryFaces(), -1);
     labelList& newFaceLabels = *newFaceLabelsPtr_;
 
-    newFaceLabelsMapPtr_ = new labelList
-    (
-        mesh_().nFaces() - mesh_().nInternalFaces(),
-        -1
-    );
+    newFaceLabelsMapPtr_ = new labelList(mesh_().nBoundaryFaces(), -1);
     labelList& newFaceLabelsMap = *newFaceLabelsMapPtr_;
     label nNewFaces = 0;
 
diff --git a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
index deb3b89bd1..f722a5578e 100644
--- a/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
+++ b/src/finiteVolume/fvMatrices/solvers/isoAdvection/isoAdvection/isoAdvection.C
@@ -92,7 +92,7 @@ Foam::isoAdvection::isoAdvection
     isoCutFace_(mesh_, ap_),
     cellIsBounded_(mesh_.nCells(), false),
     checkBounding_(mesh_.nCells(), false),
-    bsFaces_(label(0.2*(mesh_.nFaces() - mesh_.nInternalFaces()))),
+    bsFaces_(label(0.2*mesh_.nBoundaryFaces())),
     bsx0_(bsFaces_.size()),
     bsn0_(bsFaces_.size()),
     bsUn0_(bsFaces_.size()),
diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/globalIndexStencils/CFCCellToCellStencil.C b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/globalIndexStencils/CFCCellToCellStencil.C
index f7ed9ce24b..853842ee9c 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/globalIndexStencils/CFCCellToCellStencil.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/globalIndexStencils/CFCCellToCellStencil.C
@@ -36,7 +36,7 @@ void Foam::CFCCellToCellStencil::calcFaceBoundaryData
 ) const
 {
     const polyBoundaryMesh& patches = mesh().boundaryMesh();
-    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const label nBnd = mesh().nBoundaryFaces();
     const labelList& own = mesh().faceOwner();
 
     neiGlobal.setSize(nBnd);
@@ -86,7 +86,7 @@ void Foam::CFCCellToCellStencil::calcCellStencil
     labelListList& globalCellCells
 ) const
 {
-    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const label nBnd = mesh().nBoundaryFaces();
     const labelList& own = mesh().faceOwner();
     const labelList& nei = mesh().faceNeighbour();
 
diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/globalIndexStencils/cellToCellStencil.C b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/globalIndexStencils/cellToCellStencil.C
index d4b051c845..e79969b6f0 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/cellToCell/globalIndexStencils/cellToCellStencil.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/cellToCell/globalIndexStencils/cellToCellStencil.C
@@ -169,7 +169,7 @@ void Foam::cellToCellStencil::validBoundaryFaces(boolList& isValidBFace) const
 {
     const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
-    isValidBFace.setSize(mesh().nFaces()-mesh().nInternalFaces(), true);
+    isValidBFace.setSize(mesh().nBoundaryFaces(), true);
 
     for (const polyPatch& pp : patches)
     {
@@ -328,7 +328,7 @@ Foam::labelList Foam::cellToCellStencil::calcFaceCells
 Foam::cellToCellStencil::cellToCellStencil(const polyMesh& mesh)
 :
     mesh_(mesh),
-    globalNumbering_(mesh_.nCells()+mesh_.nFaces()-mesh_.nInternalFaces())
+    globalNumbering_(mesh_.nCells()+mesh_.nBoundaryFaces())
 {}
 
 
diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.C
index edd4da2e1e..eb724d3fe3 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/extendedUpwindCellToFaceStencil.C
@@ -205,7 +205,7 @@ void Foam::extendedUpwindCellToFaceStencil::transportStencils
 )
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
+    const label nBnd = mesh_.nBoundaryFaces();
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
 
diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/globalIndexStencils/FECCellToFaceStencil.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/globalIndexStencils/FECCellToFaceStencil.C
index 8012485244..ea303e3047 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/globalIndexStencils/FECCellToFaceStencil.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/globalIndexStencils/FECCellToFaceStencil.C
@@ -65,7 +65,7 @@ void Foam::FECCellToFaceStencil::calcFaceStencil
 ) const
 {
     const polyBoundaryMesh& patches = mesh().boundaryMesh();
-    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const label nBnd = mesh().nBoundaryFaces();
     const labelList& own = mesh().faceOwner();
     const labelList& nei = mesh().faceNeighbour();
 
diff --git a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/globalIndexStencils/cellToFaceStencil.C b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/globalIndexStencils/cellToFaceStencil.C
index 7106ebfae5..270047bfaa 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/cellToFace/globalIndexStencils/cellToFaceStencil.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/cellToFace/globalIndexStencils/cellToFaceStencil.C
@@ -164,7 +164,7 @@ void Foam::cellToFaceStencil::validBoundaryFaces(boolList& isValidBFace) const
 {
     const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
-    isValidBFace.setSize(mesh().nFaces()-mesh().nInternalFaces(), true);
+    isValidBFace.setSize(mesh().nBoundaryFaces(), true);
 
     forAll(patches, patchi)
     {
@@ -333,7 +333,7 @@ void Foam::cellToFaceStencil::calcFaceStencil
     // Calculates per face a list of global cell/face indices.
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-    const label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
+    const label nBnd = mesh_.nBoundaryFaces();
     const labelList& own = mesh_.faceOwner();
     const labelList& nei = mesh_.faceNeighbour();
 
@@ -482,7 +482,7 @@ void Foam::cellToFaceStencil::calcFaceStencil
 Foam::cellToFaceStencil::cellToFaceStencil(const polyMesh& mesh)
 :
     mesh_(mesh),
-    globalNumbering_(mesh_.nCells()+mesh_.nFaces()-mesh_.nInternalFaces())
+    globalNumbering_(mesh_.nCells()+mesh_.nBoundaryFaces())
 {}
 
 
diff --git a/src/finiteVolume/fvMesh/extendedStencil/faceToCell/globalIndexStencils/CFCFaceToCellStencil.C b/src/finiteVolume/fvMesh/extendedStencil/faceToCell/globalIndexStencils/CFCFaceToCellStencil.C
index 78b2b90247..efdaaaac14 100644
--- a/src/finiteVolume/fvMesh/extendedStencil/faceToCell/globalIndexStencils/CFCFaceToCellStencil.C
+++ b/src/finiteVolume/fvMesh/extendedStencil/faceToCell/globalIndexStencils/CFCFaceToCellStencil.C
@@ -36,7 +36,7 @@ void Foam::CFCFaceToCellStencil::calcFaceBoundaryData
 ) const
 {
     const polyBoundaryMesh& patches = mesh().boundaryMesh();
-    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const label nBnd = mesh().nBoundaryFaces();
     const labelList& own = mesh().faceOwner();
 
     neiGlobal.setSize(nBnd);
@@ -93,7 +93,7 @@ void Foam::CFCFaceToCellStencil::calcCellStencil
     labelListList& globalCellFaces
 ) const
 {
-    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const label nBnd = mesh().nBoundaryFaces();
     const labelList& own = mesh().faceOwner();
     const labelList& nei = mesh().faceNeighbour();
 
@@ -107,7 +107,7 @@ void Foam::CFCFaceToCellStencil::calcCellStencil
 
 
     // Non-empty boundary faces
-    boolList validBFace(mesh().nFaces()-mesh().nInternalFaces(), true);
+    boolList validBFace(mesh().nBoundaryFaces(), true);
 
     const polyBoundaryMesh& patches = mesh().boundaryMesh();
     forAll(patches, patchi)
diff --git a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C
index 1fadf245ef..d8d5b72fa5 100644
--- a/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C
+++ b/src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C
@@ -68,7 +68,7 @@ void Foam::singleCellFvMesh::agglomerateMesh
     // Check agglomeration is sync
     {
         // Get neighbouring agglomeration
-        labelList nbrAgglom(mesh.nFaces()-mesh.nInternalFaces());
+        labelList nbrAgglom(mesh.nBoundaryFaces());
         forAll(oldPatches, patchi)
         {
             const polyPatch& pp = oldPatches[patchi];
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
index a278005fda..3a2c6e12a6 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
@@ -231,7 +231,7 @@ Foam::tmp<Foam::Field<Type>> Foam::volPointInterpolation::flatBoundaryField
 
     tmp<Field<Type>> tboundaryVals
     (
-        new Field<Type>(mesh.nFaces()-mesh.nInternalFaces())
+        new Field<Type>(mesh.nBoundaryFaces())
     );
     Field<Type>& boundaryVals = tboundaryVals.ref();
 
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C
index 7b0ea37758..a36399f075 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C
@@ -57,7 +57,7 @@ void Foam::volPointInterpolation::calcBoundaryAddressing()
             SubList<face>
             (
                 mesh().faces(),
-                mesh().nFaces()-mesh().nInternalFaces(),
+                mesh().nBoundaryFaces(),
                 mesh().nInternalFaces()
             ),
             mesh().points()
diff --git a/src/functionObjects/utilities/ensightWrite/ensightWrite.C b/src/functionObjects/utilities/ensightWrite/ensightWrite.C
index 518c74c17e..75e741c2cb 100644
--- a/src/functionObjects/utilities/ensightWrite/ensightWrite.C
+++ b/src/functionObjects/utilities/ensightWrite/ensightWrite.C
@@ -214,7 +214,7 @@ bool Foam::functionObjects::ensightWrite::write()
         fileName ensightDir = dirName_;
         if (!ensightDir.isAbsolute())
         {
-            ensightDir = t.rootPath()/t.globalCaseName()/ensightDir;
+            ensightDir = t.globalPath()/ensightDir;
         }
 
         ensCase_.reset
diff --git a/src/lagrangian/basic/InteractionLists/InteractionLists.C b/src/lagrangian/basic/InteractionLists/InteractionLists.C
index 922ee09185..360f2a3cbc 100644
--- a/src/lagrangian/basic/InteractionLists/InteractionLists.C
+++ b/src/lagrangian/basic/InteractionLists/InteractionLists.C
@@ -324,10 +324,15 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
         if (isA<wallPolyPatch>(patch))
         {
-            localWallFaces.append
-            (
-                identity(patch.size(), patch.start())
-            );
+            const scalarField areaFraction(patch.areaFraction());
+
+            forAll(areaFraction, facei)
+            {
+                if (areaFraction[facei] > 0.5)
+                {
+                    localWallFaces.append(facei + patch.start());
+                }
+            }
         }
     }
 
@@ -335,10 +340,12 @@ void Foam::InteractionLists<ParticleType>::buildInteractionLists()
 
     forAll(wallFaceBbs, i)
     {
-        wallFaceBbs[i] = treeBoundBox
-        (
-            mesh_.faces()[localWallFaces[i]].points(mesh_.points())
-        );
+        wallFaceBbs[i] =
+            treeBoundBox
+            (
+                mesh_.points(),
+                mesh_.faces()[localWallFaces[i]]
+            );
     }
 
     // IAndT: index and transform
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
index da0ea0bee7..a54f1e4fa2 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.C
@@ -130,7 +130,7 @@ void Foam::meshRefinement::calcNeighbourData
     const labelList& cellLevel = meshCutter_.cellLevel();
     const pointField& cellCentres = mesh_.cellCentres();
 
-    label nBoundaryFaces = mesh_.nFaces() - mesh_.nInternalFaces();
+    const label nBoundaryFaces = mesh_.nBoundaryFaces();
 
     if (neiLevel.size() != nBoundaryFaces || neiCc.size() != nBoundaryFaces)
     {
@@ -297,8 +297,8 @@ void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
 
 
     // Get boundary face centre and level. Coupled aware.
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
+    pointField neiCc(mesh_.nBoundaryFaces());
     calcNeighbourData(neiLevel, neiCc);
 
     // Collect segments we want to test for
@@ -458,7 +458,7 @@ void Foam::meshRefinement::checkData()
     meshCutter_.checkRefinementLevels(1, labelList(0));
 
 
-    label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
+    const label nBnd = mesh_.nBoundaryFaces();
 
     Pout<< "meshRefinement::checkData() : Checking synchronization."
         << endl;
@@ -597,7 +597,7 @@ void Foam::meshRefinement::checkData()
         labelList::subList boundarySurface
         (
             surfaceIndex_,
-            mesh_.nFaces()-mesh_.nInternalFaces(),
+            mesh_.nBoundaryFaces(),
             mesh_.nInternalFaces()
         );
 
@@ -628,11 +628,7 @@ void Foam::meshRefinement::checkData()
         localPointRegion::findDuplicateFaces
         (
             mesh_,
-            identity
-            (
-                mesh_.nFaces() - mesh_.nInternalFaces(),
-                mesh_.nInternalFaces()
-            )
+            identity(mesh_.nBoundaryFaces(), mesh_.nInternalFaces())
         )
     );
 
@@ -998,12 +994,7 @@ Foam::label Foam::meshRefinement::splitFacesUndo
         // Check mesh for errors
         // ~~~~~~~~~~~~~~~~~~~~~
 
-        faceSet errorFaces
-        (
-            mesh_,
-            "errorFaces",
-            mesh_.nFaces()-mesh_.nInternalFaces()
-        );
+        faceSet errorFaces(mesh_, "errorFaces", mesh_.nBoundaryFaces());
         bool hasErrors = motionSmoother::checkMesh
         (
             false,  // report
@@ -1405,7 +1396,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::meshRefinement::balance
 
             if (keepBaffles)
             {
-                label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
+                const label nBnd = mesh_.nBoundaryFaces();
 
                 labelList coupledFace(mesh_.nFaces(), -1);
                 {
@@ -1752,7 +1743,7 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
 
     // Check that coupled faces are present on both sides.
 
-    labelList faceToZone(mesh.nFaces()-mesh.nInternalFaces(), -1);
+    labelList faceToZone(mesh.nBoundaryFaces(), -1);
 
     forAll(fZones, zonei)
     {
@@ -2891,8 +2882,8 @@ void Foam::meshRefinement::dumpIntersections(const fileName& prefix) const
         // ~~~~~~~~~~~~~~~~~~~~~~
 
         // Get boundary face centre and level. Coupled aware.
-        labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-        pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+        labelList neiLevel(mesh_.nBoundaryFaces());
+        pointField neiCc(mesh_.nBoundaryFaces());
         calcNeighbourData(neiLevel, neiCc);
 
         labelList intersectionFaces(intersectedFaces());
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
index 5ca36811c7..2fe2cbdd7a 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -2413,8 +2413,8 @@ void Foam::meshRefinement::zonify
     const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
 
     // Swap neighbouring cell centres and cell level
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
+    pointField neiCc(mesh_.nBoundaryFaces());
     calcNeighbourData(neiLevel, neiCc);
 
 
@@ -3268,7 +3268,7 @@ void Foam::meshRefinement::consistentOrientation
     {
         labelList neiStatus
         (
-            mesh_.nFaces()-mesh_.nInternalFaces(),
+            mesh_.nBoundaryFaces(),
             orientedSurface::UNVISITED
         );
 
@@ -3575,8 +3575,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
         << " faces that are intersected by the surface." << nl << endl;
 
     // Swap neighbouring cell centres and cell level
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
+    pointField neiCc(mesh_.nBoundaryFaces());
     calcNeighbourData(neiLevel, neiCc);
 
     labelList ownPatch, neiPatch;
@@ -3646,8 +3646,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
         // so re-do the surface intersections
         {
             // Swap neighbouring cell centres and cell level
-            neiLevel.setSize(mesh_.nFaces()-mesh_.nInternalFaces());
-            neiCc.setSize(mesh_.nFaces()-mesh_.nInternalFaces());
+            neiLevel.setSize(mesh_.nBoundaryFaces());
+            neiCc.setSize(mesh_.nBoundaryFaces());
             calcNeighbourData(neiLevel, neiCc);
 
             labelList ownPatch, neiPatch;
@@ -3834,8 +3834,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Swap neighbouring cell centres and cell level
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
+    pointField neiCc(mesh_.nBoundaryFaces());
     calcNeighbourData(neiLevel, neiCc);
 
     // Find intersections with all unnamed(!) surfaces
@@ -4133,8 +4133,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::removeLimitShells
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Swap neighbouring cell centres and cell level
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
+    pointField neiCc(mesh_.nBoundaryFaces());
     calcNeighbourData(neiLevel, neiCc);
 
     // Find intersections with all unnamed(!) surfaces
@@ -4459,8 +4459,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
 
     // Swap neighbouring cell centres and cell level
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
+    pointField neiCc(mesh_.nBoundaryFaces());
     calcNeighbourData(neiLevel, neiCc);
 
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
index cda200b52b..5e76cd63cc 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
@@ -90,7 +90,7 @@ Foam::label Foam::meshRefinement::markSurfaceGapRefinement
 
         // Collect cells to test for inside/outside in shell
         labelList cellToCompact(mesh_.nCells(), -1);
-        labelList bFaceToCompact(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
+        labelList bFaceToCompact(mesh_.nBoundaryFaces(), -1);
         List<FixedList<label, 3>> shellGapInfo;
         List<volumeType> shellGapMode;
         {
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
index 739bff247f..a7cc9e377e 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
@@ -49,7 +49,7 @@ Foam::label Foam::meshRefinement::mergePatchFaces
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
     // Pick up all candidate cells on boundary
-    labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelHashSet boundaryCells(mesh_.nBoundaryFaces());
 
     for (const label patchi : patchIDs)
     {
@@ -254,16 +254,14 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
     combineFaces faceCombiner(mesh_, true);
 
     // Pick up all candidate cells on boundary
-    labelHashSet boundaryCells(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelHashSet boundaryCells(mesh_.nBoundaryFaces());
 
     {
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-        forAll(patchIDs, i)
+        for (const label patchi : patchIDs)
         {
-            label patchI = patchIDs[i];
-
-            const polyPatch& patch = patches[patchI];
+            const polyPatch& patch = patches[patchi];
 
             if (!patch.coupled())
             {
@@ -410,12 +408,7 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
             // Check mesh for errors
             // ~~~~~~~~~~~~~~~~~~~~~
 
-            faceSet errorFaces
-            (
-                mesh_,
-                "errorFaces",
-                mesh_.nFaces()-mesh_.nInternalFaces()
-            );
+            faceSet errorFaces(mesh_, "errorFaces", mesh_.nBoundaryFaces());
             bool hasErrors = motionSmoother::checkMesh
             (
                 false,  // report
@@ -899,12 +892,7 @@ Foam::label Foam::meshRefinement::mergeEdgesUndo
             // Check mesh for errors
             // ~~~~~~~~~~~~~~~~~~~~~
 
-            faceSet errorFaces
-            (
-                mesh_,
-                "errorFaces",
-                mesh_.nFaces()-mesh_.nInternalFaces()
-            );
+            faceSet errorFaces(mesh_, "errorFaces", mesh_.nBoundaryFaces());
             bool hasErrors = motionSmoother::checkMesh
             (
                 false,  // report
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
index 8b147255b8..1bd63cc13c 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -520,8 +520,8 @@ Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
     labelList facePatch(mesh_.nFaces(), -1);
 
     // Swap neighbouring cell centres and cell level
-    labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiLevel(mesh_.nBoundaryFaces());
+    pointField neiCc(mesh_.nBoundaryFaces());
     calcNeighbourData(neiLevel, neiCc);
 
 
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
index 7643ab1822..182be2ee22 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
@@ -176,7 +176,7 @@ Foam::labelList Foam::meshRefinement::getChangedFaces
 
         // 2. Boundary faces
 
-        boolList refinedBoundaryFace(mesh.nFaces()-nInternalFaces, false);
+        boolList refinedBoundaryFace(mesh.nBoundaryFaces(), false);
 
         forAll(mesh.boundaryMesh(), patchI)
         {
@@ -1871,9 +1871,9 @@ Foam::label Foam::meshRefinement::markProximityRefinement
     // 2. Find out a measure of surface curvature and region edges.
     // Send over surface region and surface normal to neighbour cell.
 
-    labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
-    pointField neiBndMaxLocation(mesh_.nFaces()-mesh_.nInternalFaces());
-    vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList neiBndMaxLevel(mesh_.nBoundaryFaces());
+    pointField neiBndMaxLocation(mesh_.nBoundaryFaces());
+    vectorField neiBndMaxNormal(mesh_.nBoundaryFaces());
 
     for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
     {
diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C
index f28db75b4a..9cba909231 100644
--- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C
+++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementTemplates.C
@@ -106,7 +106,7 @@ void Foam::meshRefinement::testSyncBoundaryFaceList
     const UList<T>& syncedFaceData
 ) const
 {
-    label nBFaces = mesh_.nFaces() - mesh_.nInternalFaces();
+    const label nBFaces = mesh_.nBoundaryFaces();
 
     if (faceData.size() != nBFaces || syncedFaceData.size() != nBFaces)
     {
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
index 53b66abd60..87bb6389d4 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
@@ -165,12 +165,10 @@ void Foam::snappyLayerDriver::checkMeshManifold() const
     Info<< nl << "Checking mesh manifoldness ..." << endl;
 
     // Get all outside faces
-    labelList outsideFaces(mesh.nFaces() - mesh.nInternalFaces());
-
-    for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
-    {
-         outsideFaces[facei - mesh.nInternalFaces()] = facei;
-    }
+    labelList outsideFaces
+    (
+        identity(mesh.nBoundaryFaces(), mesh.nInternalFaces())
+    );
 
     pointSet nonManifoldPoints
     (
diff --git a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
index 34c7680aaa..a3e614313d 100644
--- a/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
+++ b/src/meshTools/AMIInterpolation/patches/cyclicPeriodicAMI/cyclicPeriodicAMIPolyPatch/cyclicPeriodicAMIPolyPatch.C
@@ -278,7 +278,7 @@ void Foam::cyclicPeriodicAMIPolyPatch::resetAMI
         {
             const Time& runTime = boundaryMesh().mesh().time();
 
-            fileName dir(runTime.rootPath()/runTime.globalCaseName());
+            fileName dir(runTime.globalPath());
             fileName postfix("_" + runTime.timeName()+"_expanded.obj");
 
             ownStr.reset(new OBJstream(dir/name() + postfix));
diff --git a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H
index 9a42db9dfa..0fbdbd5e63 100644
--- a/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H
+++ b/src/meshTools/mappedPatches/mappedPolyPatch/mappedPatchBaseI.H
@@ -111,7 +111,7 @@ inline Foam::label Foam::mappedPatchBase::sampleSize() const
         case NEARESTFACE:
         {
             const polyMesh& mesh = sampleMesh();
-            return mesh.nFaces() - mesh.nInternalFaces();
+            return mesh.nBoundaryFaces();
         }
         default:
         {
diff --git a/src/meshTools/meshSearch/meshSearch.C b/src/meshTools/meshSearch/meshSearch.C
index 9a9c7a2fd1..a14b3e841a 100644
--- a/src/meshTools/meshSearch/meshSearch.C
+++ b/src/meshTools/meshSearch/meshSearch.C
@@ -566,12 +566,11 @@ Foam::meshSearch::boundaryTree() const
             overallBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
         }
 
-        // all boundary faces (not just walls)
-        labelList bndFaces(mesh_.nFaces()-mesh_.nInternalFaces());
-        forAll(bndFaces, i)
-        {
-            bndFaces[i] = mesh_.nInternalFaces() + i;
-        }
+        // All boundary faces (not just walls)
+        labelList bndFaces
+        (
+            identity(mesh_.nBoundaryFaces(), mesh_.nInternalFaces())
+        );
 
         boundaryTreePtr_.reset
         (
diff --git a/src/meshTools/regionSplit/localPointRegion.C b/src/meshTools/regionSplit/localPointRegion.C
index 7672813d0c..5f1e0da9b2 100644
--- a/src/meshTools/regionSplit/localPointRegion.C
+++ b/src/meshTools/regionSplit/localPointRegion.C
@@ -251,7 +251,7 @@ void Foam::localPointRegion::calcPointRegions
     boolList& candidatePoint
 )
 {
-    label nBnd = mesh.nFaces()-mesh.nInternalFaces();
+    const label nBnd = mesh.nBoundaryFaces();
     const labelList& faceOwner = mesh.faceOwner();
     const labelList& faceNeighbour = mesh.faceNeighbour();
 
@@ -419,7 +419,7 @@ void Foam::localPointRegion::calcPointRegions
         SubList<face> l
         (
             minRegion,
-            mesh.nFaces()-mesh.nInternalFaces(),
+            mesh.nBoundaryFaces(),
             mesh.nInternalFaces()
         );
         syncTools::syncBoundaryFaceList
@@ -628,7 +628,7 @@ Foam::List<Foam::labelPair> Foam::localPointRegion::findDuplicateFacePairs
     // Faces to test: all boundary faces
     labelList testFaces
     (
-        identity(mesh.nFaces()-mesh.nInternalFaces(), mesh.nInternalFaces())
+        identity(mesh.nBoundaryFaces(), mesh.nInternalFaces())
     );
 
     // Find corresponding baffle face (or -1)
diff --git a/src/meshTools/regionSplit/regionSplit.C b/src/meshTools/regionSplit/regionSplit.C
index b89a64001b..47056b2092 100644
--- a/src/meshTools/regionSplit/regionSplit.C
+++ b/src/meshTools/regionSplit/regionSplit.C
@@ -530,7 +530,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::regionSplit::reduceRegions
     const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
     // Buffer for swapping boundary information
-    labelList nbrRegion(mesh().nFaces()-mesh().nInternalFaces());
+    labelList nbrRegion(mesh().nBoundaryFaces());
 
     do
     {
diff --git a/src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C b/src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C
index 7aca426d57..e4c55b966b 100644
--- a/src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C
+++ b/src/meshTools/sets/cellSources/nbrToCell/nbrToCell.C
@@ -52,7 +52,7 @@ void Foam::nbrToCell::combine(topoSet& set, const bool add) const
     const cellList& cells = mesh().cells();
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    boolList isCoupled(mesh_.nFaces()-mesh_.nInternalFaces(), false);
+    boolList isCoupled(mesh_.nBoundaryFaces(), false);
 
     for (const polyPatch& pp : patches)
     {
diff --git a/src/meshTools/sets/topoSets/faceZoneSet.C b/src/meshTools/sets/topoSets/faceZoneSet.C
index e0716e858a..674beddc35 100644
--- a/src/meshTools/sets/topoSets/faceZoneSet.C
+++ b/src/meshTools/sets/topoSets/faceZoneSet.C
@@ -351,7 +351,7 @@ void Foam::faceZoneSet::sync(const polyMesh& mesh)
     //-1 : in faceZone and flipped
     const label UNFLIPPED = 1;
     const label FLIPPED = -1;
-    labelList myZoneFace(mesh.nFaces()-mesh.nInternalFaces(), 0);
+    labelList myZoneFace(mesh.nBoundaryFaces(), Zero);
 
     forAll(addressing_, i)
     {
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 46c36f7deb..2ac294e46a 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -2283,10 +2283,7 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
     const polyMesh& mesh = bMesh.mesh();
 
     // Storage for surfaceMesh. Size estimate.
-    DynamicList<labelledTri> triangles
-    (
-        mesh.nFaces() - mesh.nInternalFaces()
-    );
+    DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
 
     label newPatchi = 0;
 
@@ -2364,10 +2361,7 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
     const polyMesh& mesh = bMesh.mesh();
 
     // Storage for surfaceMesh. Size estimate.
-    DynamicList<labelledTri> triangles
-    (
-        mesh.nFaces() - mesh.nInternalFaces()
-    );
+    DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
 
     label newPatchi = 0;
 
@@ -2472,10 +2466,7 @@ Foam::triSurface Foam::triSurfaceTools::triangulateFaceCentre
 
 
     // Count number of faces.
-    DynamicList<labelledTri> triangles
-    (
-        mesh.nFaces() - mesh.nInternalFaces()
-    );
+    DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
 
     label newPatchi = 0;
 
diff --git a/src/parallel/decompose/decompositionMethods/decompositionConstraints/preserveFaceZones/preserveFaceZonesConstraint.C b/src/parallel/decompose/decompositionMethods/decompositionConstraints/preserveFaceZones/preserveFaceZonesConstraint.C
index e481aef604..79c2400658 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionConstraints/preserveFaceZones/preserveFaceZonesConstraint.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionConstraints/preserveFaceZones/preserveFaceZonesConstraint.C
@@ -146,7 +146,7 @@ void Foam::decompositionConstraints::preserveFaceZonesConstraint::apply
 
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    labelList destProc(mesh.nFaces()-mesh.nInternalFaces(), labelMax);
+    labelList destProc(mesh.nBoundaryFaces(), labelMax);
 
     for (const polyPatch& pp : pbm)
     {
diff --git a/src/parallel/decompose/decompositionMethods/decompositionConstraints/preservePatches/preservePatchesConstraint.C b/src/parallel/decompose/decompositionMethods/decompositionConstraints/preservePatches/preservePatchesConstraint.C
index 8333cbfd7d..3513094848 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionConstraints/preservePatches/preservePatchesConstraint.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionConstraints/preservePatches/preservePatchesConstraint.C
@@ -146,7 +146,7 @@ void Foam::decompositionConstraints::preservePatchesConstraint::apply
 
     const polyBoundaryMesh& pbm = mesh.boundaryMesh();
 
-    labelList destProc(mesh.nFaces()-mesh.nInternalFaces(), labelMax);
+    labelList destProc(mesh.nBoundaryFaces(), labelMax);
 
     for (const polyPatch& pp : pbm)
     {
diff --git a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
index deaa397e03..73bb15f9a1 100644
--- a/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
+++ b/src/parallel/decompose/decompositionMethods/decompositionMethod/decompositionMethod.C
@@ -535,7 +535,7 @@ void Foam::decompositionMethod::calcCellCells
     // Get agglomerate owner on other side of coupled faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
+    labelList globalNeighbour(mesh.nBoundaryFaces());
 
     for (const polyPatch& pp : patches)
     {
@@ -734,7 +734,7 @@ void Foam::decompositionMethod::calcCellCells
     // Get agglomerate owner on other side of coupled faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    labelList globalNeighbour(mesh.nFaces()-mesh.nInternalFaces());
+    labelList globalNeighbour(mesh.nBoundaryFaces());
 
     for (const polyPatch& pp : patches)
     {
diff --git a/src/renumber/SloanRenumber/SloanRenumber.C b/src/renumber/SloanRenumber/SloanRenumber.C
index c967dd7956..075d81a89b 100644
--- a/src/renumber/SloanRenumber/SloanRenumber.C
+++ b/src/renumber/SloanRenumber/SloanRenumber.C
@@ -113,7 +113,7 @@ Foam::labelList Foam::SloanRenumber::renumber
     // Construct graph : faceOwner + connections across cyclics.
 
     // Determine neighbour cell
-    labelList nbr(mesh.nFaces()-mesh.nInternalFaces(), -1);
+    labelList nbr(mesh.nBoundaryFaces(), -1);
     forAll(pbm, patchi)
     {
         if (pbm[patchi].coupled() && !isA<processorPolyPatch>(pbm[patchi]))
diff --git a/src/sampling/meshToMesh/meshToMesh.C b/src/sampling/meshToMesh/meshToMesh.C
index 2c6a5e8ab6..f904480f00 100644
--- a/src/sampling/meshToMesh/meshToMesh.C
+++ b/src/sampling/meshToMesh/meshToMesh.C
@@ -492,7 +492,7 @@ void Foam::meshToMesh::calculate(const word& methodName, const bool normalise)
         patches[0] = new polyPatch
         (
             "defaultFaces",
-            newTgt.nFaces() - newTgt.nInternalFaces(),
+            newTgt.nBoundaryFaces(),
             newTgt.nInternalFaces(),
             0,
             newTgt.boundaryMesh(),
diff --git a/src/sampling/sampledSet/shortestPath/shortestPathSet.C b/src/sampling/sampledSet/shortestPath/shortestPathSet.C
index 2cf4276fdf..ba80d903ff 100644
--- a/src/sampling/sampledSet/shortestPath/shortestPathSet.C
+++ b/src/sampling/sampledSet/shortestPath/shortestPathSet.C
@@ -238,7 +238,7 @@ void Foam::shortestPathSet::genSamples(const polyMesh& mesh)
 
             // Situation 2: we're on a coupled patch and might need to
             //              switch processor/cell
-            boolList isFront(mesh.nFaces()-mesh.nInternalFaces(), false);
+            boolList isFront(mesh.nBoundaryFaces(), false);
 
             if (frontFaceI != -1)
             {
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 2a666f29cc..83db3ec325 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -116,7 +116,7 @@ Foam::sampledTriSurfaceMesh::nonCoupledboundaryTree() const
         // all non-coupled boundary faces (not just walls)
         const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
-        labelList bndFaces(mesh().nFaces()-mesh().nInternalFaces());
+        labelList bndFaces(patches.nFaces());
         label bndI = 0;
         for (const polyPatch& pp : patches)
         {
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
index b4f83b3736..7988cfcb03 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMeshTemplates.C
@@ -58,7 +58,7 @@ Foam::sampledTriSurfaceMesh::sampleOnFaces
     auto& values = tvalues.ref();
 
     const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-    const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+    const label nBnd = mesh().nBoundaryFaces();
 
     // Create flat boundary field
 
diff --git a/src/sampling/surface/isoSurface/isoSurface.C b/src/sampling/surface/isoSurface/isoSurface.C
index bc2a3b3c4c..1efeaec9c1 100644
--- a/src/sampling/surface/isoSurface/isoSurface.C
+++ b/src/sampling/surface/isoSurface/isoSurface.C
@@ -1453,7 +1453,7 @@ Foam::isoSurface::isoSurface
 
 
     // Pre-calculate patch-per-face to avoid whichPatch call.
-    labelList boundaryRegion(mesh_.nFaces()-mesh_.nInternalFaces());
+    labelList boundaryRegion(mesh_.nBoundaryFaces());
 
     forAll(patches, patchi)
     {
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTemplates.C b/src/sampling/surface/isoSurface/isoSurfaceTemplates.C
index 542bd28f24..c03a97de50 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTemplates.C
@@ -606,7 +606,7 @@ void Foam::isoSurface::generateTriPoints
 
 
     // Determine neighbouring snap status
-    boolList neiSnapped(mesh_.nFaces()-mesh_.nInternalFaces(), false);
+    boolList neiSnapped(mesh_.nBoundaryFaces(), false);
     List<Type> neiSnappedPoint(neiSnapped.size(), Type(Zero));
     for (const polyPatch& pp : patches)
     {
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurface.C b/src/sampling/surface/triSurfaceMesh/discreteSurface.C
index 79bfa6ddf0..5df8594ae7 100644
--- a/src/sampling/surface/triSurfaceMesh/discreteSurface.C
+++ b/src/sampling/surface/triSurfaceMesh/discreteSurface.C
@@ -110,7 +110,7 @@ Foam::discreteSurface::nonCoupledboundaryTree() const
         // all non-coupled boundary faces (not just walls)
         const polyBoundaryMesh& patches = mesh().boundaryMesh();
 
-        labelList bndFaces(mesh().nFaces()-mesh().nInternalFaces());
+        labelList bndFaces(patches.nFaces());
         label bndI = 0;
         for (const polyPatch& pp : patches)
         {
diff --git a/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C b/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
index 74256346f4..3fce29d212 100644
--- a/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
+++ b/src/sampling/surface/triSurfaceMesh/discreteSurfaceTemplates.C
@@ -169,7 +169,7 @@ Foam::discreteSurface::sampleOnFaces
         // Sample boundary faces
 
         const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-        const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+        const label nBnd = mesh().nBoundaryFaces();
 
         // Create flat boundary field
 
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index 2be602bf46..c241d13623 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.C
@@ -337,7 +337,7 @@ Foam::MeshedSurface<Face>::MeshedSurface
         SubList<face>
         (
             mesh.faces(),
-            mesh.nFaces() - mesh.nInternalFaces(),
+            mesh.nBoundaryFaces(),
             mesh.nInternalFaces()
         ),
         mesh.points()
diff --git a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
index 2c76eeeedf..4c0785eb11 100644
--- a/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
+++ b/src/thermophysicalModels/radiation/radiationModels/solarLoad/faceShading/faceShading.C
@@ -80,10 +80,7 @@ Foam::triSurface Foam::faceShading::triangulate
     const polyBoundaryMesh& bMesh = mesh_.boundaryMesh();
 
     // Storage for surfaceMesh. Size estimate.
-    DynamicList<labelledTri> triangles
-    (
-        mesh_.nFaces() - mesh_.nInternalFaces()
-    );
+    DynamicList<labelledTri> triangles(mesh_.nBoundaryFaces());
 
     label newPatchI = 0;
 
diff --git a/tutorials/incompressible/SRFPimpleFoam/rotor2D/0/epsilon b/tutorials/incompressible/SRFPimpleFoam/rotor2D/0/epsilon
index bc0943db35..e745051d37 100644
--- a/tutorials/incompressible/SRFPimpleFoam/rotor2D/0/epsilon
+++ b/tutorials/incompressible/SRFPimpleFoam/rotor2D/0/epsilon
@@ -24,7 +24,6 @@ boundaryField
     rotor
     {
         type            epsilonWallFunction;
-        U               Urel;
         value           uniform 3.75e-4;
     }
 
diff --git a/tutorials/incompressible/SRFSimpleFoam/mixer/0/epsilon b/tutorials/incompressible/SRFSimpleFoam/mixer/0/epsilon
index 63fc38dadb..5c7829c027 100644
--- a/tutorials/incompressible/SRFSimpleFoam/mixer/0/epsilon
+++ b/tutorials/incompressible/SRFSimpleFoam/mixer/0/epsilon
@@ -34,13 +34,11 @@ boundaryField
     innerWall
     {
         type            epsilonWallFunction;
-        U               Urel;
         value           uniform 14.855;
     }
     outerWall
     {
         type            epsilonWallFunction;
-        U               Urel;
         value           uniform 14.855;
     }
     cyclic_half0
diff --git a/tutorials/incompressible/SRFSimpleFoam/mixer/0/nut b/tutorials/incompressible/SRFSimpleFoam/mixer/0/nut
index 00a1844db8..2e05aa24b3 100644
--- a/tutorials/incompressible/SRFSimpleFoam/mixer/0/nut
+++ b/tutorials/incompressible/SRFSimpleFoam/mixer/0/nut
@@ -34,13 +34,11 @@ boundaryField
     innerWall
     {
         type            nutkWallFunction;
-        U               Urel;
         value           uniform 0;
     }
     outerWall
     {
         type            nutkWallFunction;
-        U               Urel;
         value           uniform 0;
     }
     cyclic_half0