diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
index bc4a363517..efeef35866 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
@@ -8,6 +8,7 @@ conformalVoronoiMesh/indexedCell/indexedCellEnum.C
 conformalVoronoiMesh/conformalVoronoiMesh.C
 conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
 conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+conformalVoronoiMesh/conformalVoronoiMeshZones.C
 conformalVoronoiMesh/conformalVoronoiMeshIO.C
 conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 9a22055dc9..eb75909d1c 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -29,10 +29,10 @@ Description
 SourceFiles
     conformalVoronoiMeshI.H
     conformalVoronoiMesh.C
+    conformalVoronoiMeshZones.C
     conformalVoronoiMeshIO.C
     conformalVoronoiMeshConformToSurface.C
     conformalVoronoiMeshFeaturePoints.C
-    conformalVoronoiMeshFeaturePointSpecialisations.C
     conformalVoronoiMeshCalcDualMesh.C
     conformalVoronoiMeshTemplates.C
 
@@ -636,6 +636,9 @@ private:
             boolList& flipMap
         ) const;
 
+        //- Add zones to the polyMesh
+        void addZones(polyMesh& mesh, const pointField& cellCentres) const;
+
         //- Tet mesh calculation
         void calcTetMesh
         (
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index dc272a89d6..4f02fb67f6 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -35,12 +35,7 @@ License
 #include "pointMesh.H"
 #include "indexedVertexOps.H"
 #include "DelaunayMeshTools.H"
-#include "surfaceZonesInfo.H"
-#include "polyModifyCell.H"
-#include "polyModifyFace.H"
 #include "syncTools.H"
-#include "regionSplit.H"
-#include "OBJstream.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -408,543 +403,6 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
 }
 
 
-void Foam::conformalVoronoiMesh::calcNeighbourCellCentres
-(
-    const polyMesh& mesh,
-    const pointField& cellCentres,
-    pointField& neiCc
-) const
-{
-    label nBoundaryFaces = mesh.nFaces() - mesh.nInternalFaces();
-
-    if (neiCc.size() != nBoundaryFaces)
-    {
-        FatalErrorIn("conformalVoronoiMesh::calcNeighbourCellCentres(..)")
-            << "nBoundaries:" << nBoundaryFaces
-            << " neiCc:" << neiCc.size()
-            << abort(FatalError);
-    }
-
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        const labelUList& faceCells = pp.faceCells();
-
-        label bFaceI = pp.start() - mesh.nInternalFaces();
-
-        if (pp.coupled())
-        {
-            forAll(faceCells, i)
-            {
-                neiCc[bFaceI] = cellCentres[faceCells[i]];
-                bFaceI++;
-            }
-        }
-    }
-
-    // Swap coupled boundaries. Apply separation to cc since is coordinate.
-    syncTools::swapBoundaryFacePositions(mesh, neiCc);
-}
-
-
-void Foam::conformalVoronoiMesh::selectSeparatedCoupledFaces
-(
-    const polyMesh& mesh,
-    boolList& selected
-) const
-{
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    forAll(patches, patchI)
-    {
-        // Check all coupled. Avoid using .coupled() so we also pick up AMI.
-        if (isA<coupledPolyPatch>(patches[patchI]))
-        {
-            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
-            (
-                patches[patchI]
-            );
-
-            if (cpp.separated() || !cpp.parallel())
-            {
-                forAll(cpp, i)
-                {
-                    selected[cpp.start()+i] = true;
-                }
-            }
-        }
-    }
-}
-
-
-void Foam::conformalVoronoiMesh::findCellZoneInsideWalk
-(
-    const polyMesh& mesh,
-    const labelList& locationSurfaces,  // indices of surfaces with inside point
-    const labelList& faceToSurface, // per face index of named surface
-    labelList& cellToSurface
-) const
-{
-    // Analyse regions. Reuse regionsplit
-    boolList blockedFace(mesh.nFaces());
-    selectSeparatedCoupledFaces(mesh, blockedFace);
-
-    forAll(faceToSurface, faceI)
-    {
-        if (faceToSurface[faceI] == -1)
-        {
-            blockedFace[faceI] = false;
-        }
-        else
-        {
-            blockedFace[faceI] = true;
-        }
-    }
-    // No need to sync since namedSurfaceIndex already is synced
-
-    // Set region per cell based on walking
-    regionSplit cellRegion(mesh, blockedFace);
-    blockedFace.clear();
-
-
-    // Force calculation of face decomposition (used in findCell)
-    (void)mesh.tetBasePtIs();
-
-    const PtrList<surfaceZonesInfo>& surfZones =
-        geometryToConformTo().surfZones();
-
-    // For all locationSurface find the cell
-    forAll(locationSurfaces, i)
-    {
-        label surfI = locationSurfaces[i];
-
-        const Foam::point& insidePoint = surfZones[surfI].zoneInsidePoint();
-
-        const word& surfName = geometryToConformTo().geometry().names()[surfI];
-
-        Info<< "    For surface " << surfName
-            << " finding inside point " << insidePoint
-            << endl;
-
-        // Find the region containing the insidePoint
-        label keepRegionI = -1;
-
-        label cellI = mesh.findCell(insidePoint);
-
-        if (cellI != -1)
-        {
-            keepRegionI = cellRegion[cellI];
-        }
-        reduce(keepRegionI, maxOp<label>());
-
-        Info<< "    For surface " << surfName
-            << " found point " << insidePoint << " in cell " << cellI
-            << " in global region " << keepRegionI
-            << " out of " << cellRegion.nRegions() << " regions." << endl;
-
-        if (keepRegionI == -1)
-        {
-            FatalErrorIn
-            (
-                "conformalVoronoiMesh::findCellZoneInsideWalk"
-                "(const polyMesh&, const labelList&"
-                ", const labelList&, labelList&)"
-            )   << "Point " << insidePoint
-                << " is not inside the mesh." << nl
-                << "Bounding box of the mesh:" << mesh.bounds()
-                << exit(FatalError);
-        }
-
-        // Set all cells with this region
-        forAll(cellRegion, cellI)
-        {
-            if (cellRegion[cellI] == keepRegionI)
-            {
-                if (cellToSurface[cellI] == -2)
-                {
-                    cellToSurface[cellI] = surfI;
-                }
-                else if (cellToSurface[cellI] != surfI)
-                {
-                    WarningIn
-                    (
-                        "conformalVoronoiMesh::findCellZoneInsideWalk"
-                        "(const labelList&, const labelList&"
-                        ", const labelList&, const labelList&)"
-                    )   << "Cell " << cellI
-                        << " at " << mesh.cellCentres()[cellI]
-                        << " is inside surface " << surfName
-                        << " but already marked as being in zone "
-                        << cellToSurface[cellI] << endl
-                        << "This can happen if your surfaces are not"
-                        << " (sufficiently) closed."
-                        << endl;
-                }
-            }
-        }
-    }
-}
-
-
-Foam::labelList Foam::conformalVoronoiMesh::calcCellZones
-(
-    const pointField& cellCentres
-) const
-{
-    labelList cellToSurface(cellCentres.size(), -1);
-
-    const PtrList<surfaceZonesInfo>& surfZones =
-        geometryToConformTo().surfZones();
-
-    // Get list of closed surfaces
-    labelList closedNamedSurfaces
-    (
-        surfaceZonesInfo::getAllClosedNamedSurfaces
-        (
-            surfZones,
-            geometryToConformTo().geometry(),
-            geometryToConformTo().surfaces()
-        )
-    );
-
-    forAll(closedNamedSurfaces, i)
-    {
-        label surfI = closedNamedSurfaces[i];
-
-        const searchableSurface& surface =
-            allGeometry()[geometryToConformTo().surfaces()[surfI]];
-
-        const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
-            surfZones[surfI].zoneInside();
-
-        if
-        (
-            selectionMethod != surfaceZonesInfo::INSIDE
-         && selectionMethod != surfaceZonesInfo::OUTSIDE
-         && selectionMethod != surfaceZonesInfo::INSIDEPOINT
-        )
-        {
-            FatalErrorIn("conformalVoronoiMesh::calcCellZones(..)")
-                << "Trying to use surface "
-                << surface.name()
-                << " which has non-geometric inside selection method "
-                << surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
-                << exit(FatalError);
-        }
-
-        if (surface.hasVolumeType())
-        {
-            List<volumeType> volType;
-            surface.getVolumeType(cellCentres, volType);
-
-            bool selectInside = true;
-            if (selectionMethod == surfaceZonesInfo::INSIDEPOINT)
-            {
-                List<volumeType> volTypeInsidePoint;
-                surface.getVolumeType
-                (
-                    pointField(1, surfZones[surfI].zoneInsidePoint()),
-                    volTypeInsidePoint
-                );
-
-                if (volTypeInsidePoint[0] == volumeType::OUTSIDE)
-                {
-                    selectInside = false;
-                }
-            }
-            else if (selectionMethod == surfaceZonesInfo::OUTSIDE)
-            {
-                selectInside = false;
-            }
-
-            forAll(volType, pointI)
-            {
-                if (cellToSurface[pointI] == -1)
-                {
-                    if
-                    (
-                        (
-                            volType[pointI] == volumeType::INSIDE
-                         && selectInside
-                        )
-                     || (
-                            volType[pointI] == volumeType::OUTSIDE
-                         && !selectInside
-                        )
-                    )
-                    {
-                        cellToSurface[pointI] = surfI;
-                    }
-                }
-            }
-        }
-    }
-
-    return cellToSurface;
-}
-
-
-void Foam::conformalVoronoiMesh::calcFaceZones
-(
-    const polyMesh& mesh,
-    const pointField& cellCentres,
-    const labelList& cellToSurface,
-    labelList& faceToSurface,
-    boolList& flipMap
-) const
-{
-    faceToSurface.setSize(mesh.nFaces(), -1);
-    flipMap.setSize(mesh.nFaces(), false);
-
-    const faceList& faces = mesh.faces();
-    const labelList& faceOwner = mesh.faceOwner();
-    const labelList& faceNeighbour = mesh.faceNeighbour();
-
-    labelList neiFaceOwner(mesh.nFaces() - mesh.nInternalFaces(), -1);
-
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        const labelUList& faceCells = pp.faceCells();
-
-        label bFaceI = pp.start() - mesh.nInternalFaces();
-
-        if (pp.coupled())
-        {
-            forAll(faceCells, i)
-            {
-                neiFaceOwner[bFaceI] = cellToSurface[faceCells[i]];
-                bFaceI++;
-            }
-        }
-    }
-
-    syncTools::swapBoundaryFaceList(mesh, neiFaceOwner);
-
-    forAll(faces, faceI)
-    {
-        const label ownerSurfaceI = cellToSurface[faceOwner[faceI]];
-
-        if (faceToSurface[faceI] >= 0)
-        {
-            continue;
-        }
-
-        if (mesh.isInternalFace(faceI))
-        {
-            const label neiSurfaceI = cellToSurface[faceNeighbour[faceI]];
-
-            if
-            (
-                (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
-             && ownerSurfaceI != neiSurfaceI
-            )
-            {
-                flipMap[faceI] =
-                    (
-                        ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
-                      ? false
-                      : true
-                    );
-
-                faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
-            }
-        }
-        else
-        {
-            label patchID = mesh.boundaryMesh().whichPatch(faceI);
-
-            if (mesh.boundaryMesh()[patchID].coupled())
-            {
-                const label neiSurfaceI =
-                    neiFaceOwner[faceI - mesh.nInternalFaces()];
-
-                if
-                (
-                    (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
-                 && ownerSurfaceI != neiSurfaceI
-                )
-                {
-                    flipMap[faceI] =
-                        (
-                            ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
-                          ? false
-                          : true
-                        );
-
-                    faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
-                }
-            }
-            else
-            {
-                if (ownerSurfaceI >= 0)
-                {
-                    faceToSurface[faceI] = ownerSurfaceI;
-                }
-            }
-        }
-    }
-
-
-    const PtrList<surfaceZonesInfo>& surfZones =
-        geometryToConformTo().surfZones();
-
-    labelList unclosedSurfaces
-    (
-        surfaceZonesInfo::getUnclosedNamedSurfaces
-        (
-            surfZones,
-            geometryToConformTo().geometry(),
-            geometryToConformTo().surfaces()
-        )
-    );
-
-    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
-    calcNeighbourCellCentres
-    (
-        mesh,
-        cellCentres,
-        neiCc
-    );
-
-    OBJstream intersections(time().path()/"ints.obj");
-    OBJstream intersectionFaces(time().path()/"intFaces.obj");
-
-    // Use intersection of cellCentre connections
-    forAll(faces, faceI)
-    {
-        if (faceToSurface[faceI] >= 0)
-        {
-            continue;
-        }
-
-        label patchID = mesh.boundaryMesh().whichPatch(faceI);
-
-        const label own = faceOwner[faceI];
-
-        List<pointIndexHit> surfHit;
-        labelList hitSurface;
-
-        if (mesh.isInternalFace(faceI))
-        {
-            const label nei = faceNeighbour[faceI];
-
-            geometryToConformTo().findSurfaceAllIntersections
-            (
-                cellCentres[own],
-                cellCentres[nei],
-                surfHit,
-                hitSurface
-            );
-        }
-        else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
-        {
-            geometryToConformTo().findSurfaceAllIntersections
-            (
-                cellCentres[own],
-                neiCc[faceI - mesh.nInternalFaces()],
-                surfHit,
-                hitSurface
-            );
-
-            if (surfHit.size() == 1 && surfHit[0].hit())
-            {
-                intersections.write
-                (
-                    linePointRef
-                    (
-                        cellCentres[own],
-                        neiCc[faceI - mesh.nInternalFaces()]
-                    )
-                );
-            }
-        }
-
-        // If there are multiple intersections then do not add to
-        // a faceZone
-        if (surfHit.size() == 1 && surfHit[0].hit())
-        {
-            if (findIndex(unclosedSurfaces, hitSurface[0]) != -1)
-            {
-                vectorField norm;
-                geometryToConformTo().getNormal
-                (
-                    hitSurface[0],
-                    List<pointIndexHit>(1, surfHit[0]),
-                    norm
-                );
-
-                vector fN = faces[faceI].normal(mesh.points());
-                fN /= mag(fN) + SMALL;
-
-                if ((norm[0] & fN) < 0)
-                {
-                    flipMap[faceI] = true;
-                }
-                else
-                {
-                    flipMap[faceI] = false;
-                }
-
-                faceToSurface[faceI] = hitSurface[0];
-
-                intersectionFaces.write(faces[faceI], mesh.points());
-            }
-        }
-    }
-
-
-//    labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
-//
-//    forAll(patches, patchI)
-//    {
-//        const polyPatch& pp = patches[patchI];
-//
-//        if (pp.coupled())
-//        {
-//            forAll(pp, i)
-//            {
-//                label faceI = pp.start()+i;
-//                label ownSurface = cellToSurface[faceOwner[faceI]];
-//                neiCellSurface[faceI - mesh.nInternalFaces()] = ownSurface;
-//            }
-//        }
-//    }
-//    syncTools::swapBoundaryFaceList(mesh, neiCellSurface);
-//
-//    forAll(patches, patchI)
-//    {
-//        const polyPatch& pp = patches[patchI];
-//
-//        if (pp.coupled())
-//        {
-//            forAll(pp, i)
-//            {
-//                label faceI = pp.start()+i;
-//                label ownSurface = cellToSurface[faceOwner[faceI]];
-//                label neiSurface =
-//                    neiCellSurface[faceI-mesh.nInternalFaces()];
-//
-//                if (faceToSurface[faceI] == -1 && (ownSurface != neiSurface))
-//                {
-//                    // Give face the max cell zone
-//                    faceToSurface[faceI] =  max(ownSurface, neiSurface);
-//                }
-//            }
-//        }
-//    }
-//
-    // Sync
-    syncTools::syncFaceList(mesh, faceToSurface, maxEqOp<label>());
-}
-
-
 Foam::autoPtr<Foam::fvMesh> Foam::conformalVoronoiMesh::createDummyMesh
 (
     const IOobject& io,
@@ -1304,7 +762,7 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
         {
             if (rotation[faceI] != 0)
             {
-                inplaceRotateList<List, label>(faces[faceI], rotation[faceI]);
+                faces[faceI] = rotateList(faces[faceI], rotation[faceI]);
                 nRotated++;
             }
         }
@@ -1448,159 +906,9 @@ void Foam::conformalVoronoiMesh::writeMesh
     mesh.addFvPatches(patches);
 
 
-    // Add zones to mesh
-    {
-        Info<< "    Adding zones to mesh" << endl;
-
-        const PtrList<surfaceZonesInfo>& surfZones =
-            geometryToConformTo().surfZones();
-
-        labelList cellToSurface(calcCellZones(cellCentres));
-
-        labelList faceToSurface;
-        boolList flipMap;
-
-        calcFaceZones
-        (
-            mesh,
-            cellCentres,
-            cellToSurface,
-            faceToSurface,
-            flipMap
-        );
-
-        labelList insidePointNamedSurfaces
-        (
-            surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
-        );
-
-        findCellZoneInsideWalk
-        (
-            mesh,
-            insidePointNamedSurfaces,
-            faceToSurface,
-            cellToSurface
-        );
-
-        labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
-
-        forAll(namedSurfaces, i)
-        {
-            label surfI = namedSurfaces[i];
-
-            Info<< incrIndent << indent << "Surface : "
-                << geometryToConformTo().geometry().names()[surfI] << nl
-                << indent << "    faceZone : "
-                << surfZones[surfI].faceZoneName() << nl
-                << indent << "    cellZone : "
-                << surfZones[surfI].cellZoneName()
-                << decrIndent << endl;
-        }
-
-        // Add zones to mesh
-        labelList surfaceToFaceZone =
-            surfaceZonesInfo::addFaceZonesToMesh
-            (
-                surfZones,
-                namedSurfaces,
-                mesh
-            );
-
-        labelList surfaceToCellZone =
-            surfaceZonesInfo::addCellZonesToMesh
-            (
-                surfZones,
-                namedSurfaces,
-                mesh
-            );
-
-        // Topochange container
-        polyTopoChange meshMod(mesh);
-
-        forAll(cellToSurface, cellI)
-        {
-            label surfaceI = cellToSurface[cellI];
-
-            if (surfaceI >= 0)
-            {
-                label zoneI = surfaceToCellZone[surfaceI];
-
-                if (zoneI >= 0)
-                {
-                    meshMod.setAction
-                    (
-                        polyModifyCell
-                        (
-                            cellI,
-                            false,          // removeFromZone
-                            zoneI
-                        )
-                    );
-                }
-            }
-        }
-
-        const labelList& faceOwner = mesh.faceOwner();
-        const labelList& faceNeighbour = mesh.faceNeighbour();
-
-        forAll(faceToSurface, faceI)
-        {
-            label surfaceI = faceToSurface[faceI];
-
-            if (surfaceI < 0)
-            {
-                continue;
-            }
-
-            label patchID = mesh.boundaryMesh().whichPatch(faceI);
-
-            if (mesh.isInternalFace(faceI))
-            {
-                label own = faceOwner[faceI];
-                label nei = faceNeighbour[faceI];
-
-                meshMod.setAction
-                (
-                    polyModifyFace
-                    (
-                        mesh.faces()[faceI],            // modified face
-                        faceI,                          // label of face
-                        own,                            // owner
-                        nei,                            // neighbour
-                        false,                          // face flip
-                        -1,                             // patch for face
-                        false,                          // remove from zone
-                        surfaceToFaceZone[surfaceI],    // zone for face
-                        flipMap[faceI]                  // face flip in zone
-                    )
-                );
-            }
-            else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
-            {
-                label own = faceOwner[faceI];
-
-                meshMod.setAction
-                (
-                    polyModifyFace
-                    (
-                        mesh.faces()[faceI],            // modified face
-                        faceI,                          // label of face
-                        own,                            // owner
-                        -1,                             // neighbour
-                        false,                          // face flip
-                        patchID,                        // patch for face
-                        false,                          // remove from zone
-                        surfaceToFaceZone[surfaceI],    // zone for face
-                        flipMap[faceI]                  // face flip in zone
-                    )
-                );
-            }
-        }
-
-        // Change the mesh (no inflation, parallel sync)
-        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
-    }
 
+    // Add zones to the mesh
+    addZones(mesh, cellCentres);
 
 
 
@@ -1611,12 +919,18 @@ void Foam::conformalVoronoiMesh::writeMesh
 
         forAll(boundaryFacesToRemove, faceI)
         {
-            if (boundaryFacesToRemove[faceI])
+            if
+            (
+                boundaryFacesToRemove[faceI]
+             && mesh.faceZones().whichZone(faceI) == -1
+            )
             {
                 addr[count++] = faceI;
             }
         }
 
+        addr.setSize(count);
+
         label sz = mesh.faceZones().size();
         boolList flip(addr.size(), false);
         mesh.faceZones().setSize(sz + 1);
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C
new file mode 100644
index 0000000000..e0915e01fb
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C
@@ -0,0 +1,713 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "conformalVoronoiMesh.H"
+#include "polyModifyFace.H"
+#include "polyModifyCell.H"
+#include "syncTools.H"
+#include "regionSplit.H"
+#include "surfaceZonesInfo.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::conformalVoronoiMesh::calcNeighbourCellCentres
+(
+    const polyMesh& mesh,
+    const pointField& cellCentres,
+    pointField& neiCc
+) const
+{
+    label nBoundaryFaces = mesh.nFaces() - mesh.nInternalFaces();
+
+    if (neiCc.size() != nBoundaryFaces)
+    {
+        FatalErrorIn("conformalVoronoiMesh::calcNeighbourCellCentres(..)")
+            << "nBoundaries:" << nBoundaryFaces
+            << " neiCc:" << neiCc.size()
+            << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        const labelUList& faceCells = pp.faceCells();
+
+        label bFaceI = pp.start() - mesh.nInternalFaces();
+
+        if (pp.coupled())
+        {
+            forAll(faceCells, i)
+            {
+                neiCc[bFaceI] = cellCentres[faceCells[i]];
+                bFaceI++;
+            }
+        }
+    }
+
+    // Swap coupled boundaries. Apply separation to cc since is coordinate.
+    syncTools::swapBoundaryFacePositions(mesh, neiCc);
+}
+
+
+void Foam::conformalVoronoiMesh::selectSeparatedCoupledFaces
+(
+    const polyMesh& mesh,
+    boolList& selected
+) const
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        // Check all coupled. Avoid using .coupled() so we also pick up AMI.
+        if (isA<coupledPolyPatch>(patches[patchI]))
+        {
+            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+            (
+                patches[patchI]
+            );
+
+            if (cpp.separated() || !cpp.parallel())
+            {
+                forAll(cpp, i)
+                {
+                    selected[cpp.start()+i] = true;
+                }
+            }
+        }
+    }
+}
+
+
+void Foam::conformalVoronoiMesh::findCellZoneInsideWalk
+(
+    const polyMesh& mesh,
+    const labelList& locationSurfaces,  // indices of surfaces with inside point
+    const labelList& faceToSurface, // per face index of named surface
+    labelList& cellToSurface
+) const
+{
+    // Analyse regions. Reuse regionsplit
+    boolList blockedFace(mesh.nFaces());
+    selectSeparatedCoupledFaces(mesh, blockedFace);
+
+    forAll(faceToSurface, faceI)
+    {
+        if (faceToSurface[faceI] == -1)
+        {
+            blockedFace[faceI] = false;
+        }
+        else
+        {
+            blockedFace[faceI] = true;
+        }
+    }
+    // No need to sync since namedSurfaceIndex already is synced
+
+    // Set region per cell based on walking
+    regionSplit cellRegion(mesh, blockedFace);
+    blockedFace.clear();
+
+
+    // Force calculation of face decomposition (used in findCell)
+    (void)mesh.tetBasePtIs();
+
+    const PtrList<surfaceZonesInfo>& surfZones =
+        geometryToConformTo().surfZones();
+
+    // For all locationSurface find the cell
+    forAll(locationSurfaces, i)
+    {
+        label surfI = locationSurfaces[i];
+
+        const Foam::point& insidePoint = surfZones[surfI].zoneInsidePoint();
+
+        const word& surfName = geometryToConformTo().geometry().names()[surfI];
+
+        Info<< "    For surface " << surfName
+            << " finding inside point " << insidePoint
+            << endl;
+
+        // Find the region containing the insidePoint
+        label keepRegionI = -1;
+
+        label cellI = mesh.findCell(insidePoint);
+
+        if (cellI != -1)
+        {
+            keepRegionI = cellRegion[cellI];
+        }
+        reduce(keepRegionI, maxOp<label>());
+
+        Info<< "    For surface " << surfName
+            << " found point " << insidePoint << " in cell " << cellI
+            << " in global region " << keepRegionI
+            << " out of " << cellRegion.nRegions() << " regions." << endl;
+
+        if (keepRegionI == -1)
+        {
+            FatalErrorIn
+            (
+                "conformalVoronoiMesh::findCellZoneInsideWalk"
+                "(const polyMesh&, const labelList&"
+                ", const labelList&, labelList&)"
+            )   << "Point " << insidePoint
+                << " is not inside the mesh." << nl
+                << "Bounding box of the mesh:" << mesh.bounds()
+                << exit(FatalError);
+        }
+
+        // Set all cells with this region
+        forAll(cellRegion, cellI)
+        {
+            if (cellRegion[cellI] == keepRegionI)
+            {
+                if (cellToSurface[cellI] == -2)
+                {
+                    cellToSurface[cellI] = surfI;
+                }
+                else if (cellToSurface[cellI] != surfI)
+                {
+                    WarningIn
+                    (
+                        "conformalVoronoiMesh::findCellZoneInsideWalk"
+                        "(const labelList&, const labelList&"
+                        ", const labelList&, const labelList&)"
+                    )   << "Cell " << cellI
+                        << " at " << mesh.cellCentres()[cellI]
+                        << " is inside surface " << surfName
+                        << " but already marked as being in zone "
+                        << cellToSurface[cellI] << endl
+                        << "This can happen if your surfaces are not"
+                        << " (sufficiently) closed."
+                        << endl;
+                }
+            }
+        }
+    }
+}
+
+
+Foam::labelList Foam::conformalVoronoiMesh::calcCellZones
+(
+    const pointField& cellCentres
+) const
+{
+    labelList cellToSurface(cellCentres.size(), -1);
+
+    const PtrList<surfaceZonesInfo>& surfZones =
+        geometryToConformTo().surfZones();
+
+    // Get list of closed surfaces
+    labelList closedNamedSurfaces
+    (
+        surfaceZonesInfo::getAllClosedNamedSurfaces
+        (
+            surfZones,
+            geometryToConformTo().geometry(),
+            geometryToConformTo().surfaces()
+        )
+    );
+
+    forAll(closedNamedSurfaces, i)
+    {
+        label surfI = closedNamedSurfaces[i];
+
+        const searchableSurface& surface =
+            allGeometry()[geometryToConformTo().surfaces()[surfI]];
+
+        const surfaceZonesInfo::areaSelectionAlgo selectionMethod =
+            surfZones[surfI].zoneInside();
+
+        if
+        (
+            selectionMethod != surfaceZonesInfo::INSIDE
+         && selectionMethod != surfaceZonesInfo::OUTSIDE
+         && selectionMethod != surfaceZonesInfo::INSIDEPOINT
+        )
+        {
+            FatalErrorIn("conformalVoronoiMesh::calcCellZones(..)")
+                << "Trying to use surface "
+                << surface.name()
+                << " which has non-geometric inside selection method "
+                << surfaceZonesInfo::areaSelectionAlgoNames[selectionMethod]
+                << exit(FatalError);
+        }
+
+        if (surface.hasVolumeType())
+        {
+            List<volumeType> volType;
+            surface.getVolumeType(cellCentres, volType);
+
+            bool selectInside = true;
+            if (selectionMethod == surfaceZonesInfo::INSIDEPOINT)
+            {
+                List<volumeType> volTypeInsidePoint;
+                surface.getVolumeType
+                (
+                    pointField(1, surfZones[surfI].zoneInsidePoint()),
+                    volTypeInsidePoint
+                );
+
+                if (volTypeInsidePoint[0] == volumeType::OUTSIDE)
+                {
+                    selectInside = false;
+                }
+            }
+            else if (selectionMethod == surfaceZonesInfo::OUTSIDE)
+            {
+                selectInside = false;
+            }
+
+            forAll(volType, pointI)
+            {
+                if (cellToSurface[pointI] == -1)
+                {
+                    if
+                    (
+                        (
+                            volType[pointI] == volumeType::INSIDE
+                         && selectInside
+                        )
+                     || (
+                            volType[pointI] == volumeType::OUTSIDE
+                         && !selectInside
+                        )
+                    )
+                    {
+                        cellToSurface[pointI] = surfI;
+                    }
+                }
+            }
+        }
+    }
+
+    return cellToSurface;
+}
+
+
+void Foam::conformalVoronoiMesh::calcFaceZones
+(
+    const polyMesh& mesh,
+    const pointField& cellCentres,
+    const labelList& cellToSurface,
+    labelList& faceToSurface,
+    boolList& flipMap
+) const
+{
+    faceToSurface.setSize(mesh.nFaces(), -1);
+    flipMap.setSize(mesh.nFaces(), false);
+
+    const faceList& faces = mesh.faces();
+    const labelList& faceOwner = mesh.faceOwner();
+    const labelList& faceNeighbour = mesh.faceNeighbour();
+
+    labelList neiFaceOwner(mesh.nFaces() - mesh.nInternalFaces(), -1);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        const labelUList& faceCells = pp.faceCells();
+
+        label bFaceI = pp.start() - mesh.nInternalFaces();
+
+        if (pp.coupled())
+        {
+            forAll(faceCells, i)
+            {
+                neiFaceOwner[bFaceI] = cellToSurface[faceCells[i]];
+                bFaceI++;
+            }
+        }
+    }
+
+    syncTools::swapBoundaryFaceList(mesh, neiFaceOwner);
+
+    forAll(faces, faceI)
+    {
+        const label ownerSurfaceI = cellToSurface[faceOwner[faceI]];
+
+        if (faceToSurface[faceI] >= 0)
+        {
+            continue;
+        }
+
+        if (mesh.isInternalFace(faceI))
+        {
+            const label neiSurfaceI = cellToSurface[faceNeighbour[faceI]];
+
+            if
+            (
+                (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
+             && ownerSurfaceI != neiSurfaceI
+            )
+            {
+                flipMap[faceI] =
+                    (
+                        ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
+                      ? false
+                      : true
+                    );
+
+                faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
+            }
+        }
+        else
+        {
+            label patchID = mesh.boundaryMesh().whichPatch(faceI);
+
+            if (mesh.boundaryMesh()[patchID].coupled())
+            {
+                const label neiSurfaceI =
+                    neiFaceOwner[faceI - mesh.nInternalFaces()];
+
+                if
+                (
+                    (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
+                 && ownerSurfaceI != neiSurfaceI
+                )
+                {
+                    flipMap[faceI] =
+                        (
+                            ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
+                          ? false
+                          : true
+                        );
+
+                    faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
+                }
+            }
+            else
+            {
+                if (ownerSurfaceI >= 0)
+                {
+                    faceToSurface[faceI] = ownerSurfaceI;
+                }
+            }
+        }
+    }
+
+
+    const PtrList<surfaceZonesInfo>& surfZones =
+        geometryToConformTo().surfZones();
+
+    labelList unclosedSurfaces
+    (
+        surfaceZonesInfo::getUnclosedNamedSurfaces
+        (
+            surfZones,
+            geometryToConformTo().geometry(),
+            geometryToConformTo().surfaces()
+        )
+    );
+
+    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
+    calcNeighbourCellCentres
+    (
+        mesh,
+        cellCentres,
+        neiCc
+    );
+
+    // Use intersection of cellCentre connections
+    forAll(faces, faceI)
+    {
+        if (faceToSurface[faceI] >= 0)
+        {
+            continue;
+        }
+
+        label patchID = mesh.boundaryMesh().whichPatch(faceI);
+
+        const label own = faceOwner[faceI];
+
+        List<pointIndexHit> surfHit;
+        labelList hitSurface;
+
+        if (mesh.isInternalFace(faceI))
+        {
+            const label nei = faceNeighbour[faceI];
+
+            geometryToConformTo().findSurfaceAllIntersections
+            (
+                cellCentres[own],
+                cellCentres[nei],
+                surfHit,
+                hitSurface
+            );
+        }
+        else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
+        {
+            geometryToConformTo().findSurfaceAllIntersections
+            (
+                cellCentres[own],
+                neiCc[faceI - mesh.nInternalFaces()],
+                surfHit,
+                hitSurface
+            );
+        }
+
+        // If there are multiple intersections then do not add to
+        // a faceZone
+        if (surfHit.size() == 1 && surfHit[0].hit())
+        {
+            if (findIndex(unclosedSurfaces, hitSurface[0]) != -1)
+            {
+                vectorField norm;
+                geometryToConformTo().getNormal
+                (
+                    hitSurface[0],
+                    List<pointIndexHit>(1, surfHit[0]),
+                    norm
+                );
+
+                vector fN = faces[faceI].normal(mesh.points());
+                fN /= mag(fN) + SMALL;
+
+                if ((norm[0] & fN) < 0)
+                {
+                    flipMap[faceI] = true;
+                }
+                else
+                {
+                    flipMap[faceI] = false;
+                }
+
+                faceToSurface[faceI] = hitSurface[0];
+            }
+        }
+    }
+
+
+//    labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
+//
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (pp.coupled())
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//                label ownSurface = cellToSurface[faceOwner[faceI]];
+//                neiCellSurface[faceI - mesh.nInternalFaces()] = ownSurface;
+//            }
+//        }
+//    }
+//    syncTools::swapBoundaryFaceList(mesh, neiCellSurface);
+//
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (pp.coupled())
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//                label ownSurface = cellToSurface[faceOwner[faceI]];
+//                label neiSurface =
+//                    neiCellSurface[faceI-mesh.nInternalFaces()];
+//
+//                if (faceToSurface[faceI] == -1 && (ownSurface != neiSurface))
+//                {
+//                    // Give face the max cell zone
+//                    faceToSurface[faceI] =  max(ownSurface, neiSurface);
+//                }
+//            }
+//        }
+//    }
+
+    // Sync
+    syncTools::syncFaceList(mesh, faceToSurface, maxEqOp<label>());
+}
+
+
+void Foam::conformalVoronoiMesh::addZones
+(
+    polyMesh& mesh,
+    const pointField& cellCentres
+) const
+{
+    Info<< "    Adding zones to mesh" << endl;
+
+    const PtrList<surfaceZonesInfo>& surfZones =
+        geometryToConformTo().surfZones();
+
+    labelList cellToSurface(calcCellZones(cellCentres));
+
+    labelList faceToSurface;
+    boolList flipMap;
+
+    calcFaceZones
+    (
+        mesh,
+        cellCentres,
+        cellToSurface,
+        faceToSurface,
+        flipMap
+    );
+
+    labelList insidePointNamedSurfaces
+    (
+        surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
+    );
+
+    findCellZoneInsideWalk
+    (
+        mesh,
+        insidePointNamedSurfaces,
+        faceToSurface,
+        cellToSurface
+    );
+
+    labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
+
+    forAll(namedSurfaces, i)
+    {
+        label surfI = namedSurfaces[i];
+
+        Info<< incrIndent << indent << "Surface : "
+            << geometryToConformTo().geometry().names()[surfI] << nl
+            << indent << "    faceZone : "
+            << surfZones[surfI].faceZoneName() << nl
+            << indent << "    cellZone : "
+            << surfZones[surfI].cellZoneName()
+            << decrIndent << endl;
+    }
+
+    // Add zones to mesh
+    labelList surfaceToFaceZone =
+        surfaceZonesInfo::addFaceZonesToMesh
+        (
+            surfZones,
+            namedSurfaces,
+            mesh
+        );
+
+    labelList surfaceToCellZone =
+        surfaceZonesInfo::addCellZonesToMesh
+        (
+            surfZones,
+            namedSurfaces,
+            mesh
+        );
+
+    // Topochange container
+    polyTopoChange meshMod(mesh);
+
+    forAll(cellToSurface, cellI)
+    {
+        label surfaceI = cellToSurface[cellI];
+
+        if (surfaceI >= 0)
+        {
+            label zoneI = surfaceToCellZone[surfaceI];
+
+            if (zoneI >= 0)
+            {
+                meshMod.setAction
+                (
+                    polyModifyCell
+                    (
+                        cellI,
+                        false,          // removeFromZone
+                        zoneI
+                    )
+                );
+            }
+        }
+    }
+
+    const labelList& faceOwner = mesh.faceOwner();
+    const labelList& faceNeighbour = mesh.faceNeighbour();
+
+    forAll(faceToSurface, faceI)
+    {
+        label surfaceI = faceToSurface[faceI];
+
+        if (surfaceI < 0)
+        {
+            continue;
+        }
+
+        label patchID = mesh.boundaryMesh().whichPatch(faceI);
+
+        if (mesh.isInternalFace(faceI))
+        {
+            label own = faceOwner[faceI];
+            label nei = faceNeighbour[faceI];
+
+            meshMod.setAction
+            (
+                polyModifyFace
+                (
+                    mesh.faces()[faceI],            // modified face
+                    faceI,                          // label of face
+                    own,                            // owner
+                    nei,                            // neighbour
+                    false,                          // face flip
+                    -1,                             // patch for face
+                    false,                          // remove from zone
+                    surfaceToFaceZone[surfaceI],    // zone for face
+                    flipMap[faceI]                  // face flip in zone
+                )
+            );
+        }
+        else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
+        {
+            label own = faceOwner[faceI];
+
+            meshMod.setAction
+            (
+                polyModifyFace
+                (
+                    mesh.faces()[faceI],            // modified face
+                    faceI,                          // label of face
+                    own,                            // owner
+                    -1,                             // neighbour
+                    false,                          // face flip
+                    patchID,                        // patch for face
+                    false,                          // remove from zone
+                    surfaceToFaceZone[surfaceI],    // zone for face
+                    flipMap[faceI]                  // face flip in zone
+                )
+            );
+        }
+    }
+
+    // Change the mesh (no inflation, parallel sync)
+    autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H b/src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H
index f6e3f2a585..fe0e33da2d 100644
--- a/src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H
+++ b/src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,7 +36,7 @@ Description
 
     Storage is allocated on free-store during construction.
 
-    As a special case a null-contructed CompactListList has an empty
+    As a special case a null-constructed CompactListList has an empty
     offsets_ (instead of size 1).
 
 SourceFiles
diff --git a/src/OpenFOAM/db/IOobject/IOobject.C b/src/OpenFOAM/db/IOobject/IOobject.C
index 65fddffbbc..f19b113b55 100644
--- a/src/OpenFOAM/db/IOobject/IOobject.C
+++ b/src/OpenFOAM/db/IOobject/IOobject.C
@@ -63,23 +63,23 @@ bool Foam::IOobject::IOobject::fileNameComponents
     // called with directory
     if (isDir(path))
     {
-        WarningIn("IOobject::fileNameComponents(const fileName&, ...)")
-            << " called with directory: " << path << "\n";
+        WarningIn
+        (
+            "IOobject::fileNameComponents"
+            "("
+                "const fileName&, "
+                "fileName&, "
+                "fileName&, "
+                "word&"
+            ")"
+        )
+            << " called with directory: " << path << endl;
+
         return false;
     }
 
-    string::size_type first = path.find('/');
-
-    if (first == string::npos)
+    if (path.isAbsolute())
     {
-        // no '/' found - no instance or local
-
-        // check afterwards
-        name.string::operator=(path);
-    }
-    else if (first == 0)
-    {
-        // Leading '/'. Absolute fileName
         string::size_type last = path.rfind('/');
         instance = path.substr(0, last);
 
@@ -88,26 +88,48 @@ bool Foam::IOobject::IOobject::fileNameComponents
     }
     else
     {
-        instance = path.substr(0, first);
+        string::size_type first = path.find('/');
 
-        string::size_type last = path.rfind('/');
-        if (last > first)
+        if (first == string::npos)
         {
-            // with local
-            local = path.substr(first+1, last-first-1);
-        }
+            // no '/' found - no instance or local
 
-        // check afterwards
-        name.string::operator=(path.substr(last+1));
+            // check afterwards
+            name.string::operator=(path);
+        }
+        else
+        {
+            instance = path.substr(0, first);
+
+            string::size_type last = path.rfind('/');
+            if (last > first)
+            {
+                // with local
+                local = path.substr(first+1, last-first-1);
+            }
+
+            // check afterwards
+            name.string::operator=(path.substr(last+1));
+        }
     }
 
 
     // check for valid (and stripped) name, regardless of the debug level
     if (name.empty() || string::stripInvalid<word>(name))
     {
-        WarningIn("IOobject::fileNameComponents(const fileName&, ...)")
+        WarningIn
+        (
+            "IOobject::fileNameComponents"
+            "("
+                "const fileName&, "
+                "fileName&, "
+                "fileName&, "
+                "word&"
+            ")"
+        )
             << "has invalid word for name: \"" << name
-            << "\"\nwhile processing path: " << path << "\n";
+            << "\"\nwhile processing path: " << path << endl;
+
         return false;
     }
 
@@ -202,9 +224,16 @@ Foam::IOobject::IOobject
     {
         FatalErrorIn
         (
-            "IOobject::IOobject" "(const fileName&, const objectRegistry&, ...)"
+            "IOobject::IOobject"
+            "("
+                "const fileName&, "
+                "const objectRegistry&, "
+                "readOption, "
+                "writeOption, "
+                "bool"
+            ")"
         )
-            << " invalid path specification\n"
+            << " invalid path specification"
             << exit(FatalError);
     }
 
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
index 905a565b0a..9fa08201ea 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C
@@ -294,15 +294,26 @@ void Foam::cyclicACMIFvPatchField<Type>::updateCoeffs()
 }
 
 
+template<class Type>
+void Foam::cyclicACMIFvPatchField<Type>::initEvaluate
+(
+    const Pstream::commsTypes comms
+)
+{
+    // update non-overlap patch
+    const fvPatchField<Type>& npf = nonOverlapPatchField();
+    const_cast<fvPatchField<Type>&>(npf).evaluate(comms);
+}
+
+
 template<class Type>
 void Foam::cyclicACMIFvPatchField<Type>::evaluate
 (
     const Pstream::commsTypes comms
 )
 {
-    // blend contrubutions from the coupled and non-overlap patches
+    // blend contributions from the coupled and non-overlap patches
     const fvPatchField<Type>& npf = nonOverlapPatchField();
-    const_cast<fvPatchField<Type>&>(npf).evaluate();
 
     coupledFvPatchField<Type>::evaluate(comms);
     const Field<Type>& cpf = *this;
@@ -392,7 +403,7 @@ void Foam::cyclicACMIFvPatchField<Type>::manipulateMatrix
     fvMatrix<Type>& matrix
 )
 {
-    // blend contrubutions from the coupled and non-overlap patches
+    // blend contributions from the coupled and non-overlap patches
     const fvPatchField<Type>& npf = nonOverlapPatchField();
 
     const scalarField& mask = cyclicACMIPatch_.cyclicACMIPatch().mask();
diff --git a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
index 91f414d321..3a8198a95b 100644
--- a/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.H
@@ -206,6 +206,12 @@ public:
             //- Update the coefficients associated with the patch field
             void updateCoeffs();
 
+            //- Initialise the evaluation of the patch field
+            virtual void initEvaluate
+            (
+                const Pstream::commsTypes commsType
+            );
+
             //- Evaluate the patch field
             virtual void evaluate
             (
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
index c29dfe2c8d..e32a3bdf4f 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H
@@ -397,6 +397,9 @@ public:
 //            inline const typename parcelType::forceType& forces() const;
             inline const forceType& forces() const;
 
+            //- Return the optional particle forces
+            inline forceType& forces();
+
             //- Optional cloud function objects
             inline functionType& functions();
 
@@ -487,6 +490,9 @@ public:
             //- Total linear kinetic energy in the system
             inline scalar linearKineticEnergyOfSystem() const;
 
+            //- Total rotational kinetic energy in the system
+            inline scalar rotationalKineticEnergyOfSystem() const;
+
             //- Penetration for fraction [0-1] of the current total mass
             inline scalar penetration(const scalar fraction) const;
 
diff --git a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
index 5f0d6767f3..1bbf47b41a 100644
--- a/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
+++ b/src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H
@@ -156,6 +156,14 @@ Foam::KinematicCloud<CloudType>::forces() const
 }
 
 
+template<class CloudType>
+inline typename Foam::KinematicCloud<CloudType>::forceType&
+Foam::KinematicCloud<CloudType>::forces()
+{
+    return forces_;
+}
+
+
 template<class CloudType>
 inline typename Foam::KinematicCloud<CloudType>::functionType&
 Foam::KinematicCloud<CloudType>::functions()
diff --git a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C
index 04016823b4..81fced59a8 100644
--- a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C
+++ b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,9 @@ Foam::ParticleForceList<CloudType>::ParticleForceList
     PtrList<ParticleForce<CloudType> >(),
     owner_(owner),
     mesh_(mesh),
-    dict_(dictionary::null)
+    dict_(dictionary::null),
+    calcCoupled_(true),
+    calcNonCoupled_(true)
 {}
 
 
@@ -54,7 +56,9 @@ Foam::ParticleForceList<CloudType>::ParticleForceList
     PtrList<ParticleForce<CloudType> >(),
     owner_(owner),
     mesh_(mesh),
-    dict_(dict)
+    dict_(dict),
+    calcCoupled_(true),
+    calcNonCoupled_(true)
 {
     if (readFields)
     {
@@ -151,9 +155,13 @@ Foam::forceSuSp Foam::ParticleForceList<CloudType>::calcCoupled
 ) const
 {
     forceSuSp value(vector::zero, 0.0);
-    forAll(*this, i)
+
+    if (calcCoupled_)
     {
-        value += this->operator[](i).calcCoupled(p, dt, mass, Re, muc);
+        forAll(*this, i)
+        {
+            value += this->operator[](i).calcCoupled(p, dt, mass, Re, muc);
+        }
     }
 
     return value;
@@ -171,9 +179,13 @@ Foam::forceSuSp Foam::ParticleForceList<CloudType>::calcNonCoupled
 ) const
 {
     forceSuSp value(vector::zero, 0.0);
-    forAll(*this, i)
+
+    if (calcNonCoupled_)
     {
-        value += this->operator[](i).calcNonCoupled(p, dt, mass, Re, muc);
+        forAll(*this, i)
+        {
+            value += this->operator[](i).calcNonCoupled(p, dt, mass, Re, muc);
+        }
     }
 
     return value;
diff --git a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
index fe5f561268..0c7b08a7ed 100644
--- a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
+++ b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceList.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,6 +64,12 @@ class ParticleForceList
         //- Forces dictionary
         const dictionary dict_;
 
+        //- Calculate coupled forces flag
+        bool calcCoupled_;
+
+        //- Calculate non-coupled forces flag
+        bool calcNonCoupled_;
+
 
 public:
 
@@ -105,6 +111,12 @@ public:
             //- Return the forces dictionary
             inline const dictionary& dict() const;
 
+            //- Set the calcCoupled flag
+            inline void setCalcCoupled(bool flag);
+
+            //- Set the calcNonCoupled flag
+            inline void setCalcNonCoupled(bool flag);
+
 
         // Evaluation
 
diff --git a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceListI.H b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceListI.H
index db2f807ed2..a803255935 100644
--- a/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceListI.H
+++ b/src/lagrangian/intermediate/submodels/ForceTypes/ParticleForceList/ParticleForceListI.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,4 +53,18 @@ inline const Foam::dictionary& Foam::ParticleForceList<CloudType>::dict() const
 }
 
 
+template<class CloudType>
+inline void Foam::ParticleForceList<CloudType>::setCalcCoupled(bool flag)
+{
+    calcCoupled_ = flag;
+}
+
+
+template<class CloudType>
+inline void Foam::ParticleForceList<CloudType>::setCalcNonCoupled(bool flag)
+{
+    calcNonCoupled_ = flag;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
index 198dc717f8..563a6763a2 100644
--- a/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
+++ b/src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C
@@ -68,13 +68,11 @@ void Foam::SprayParcel<ParcelType>::calc
     const CompositionModel<reactingCloudType>& composition =
         td.cloud().composition();
 
-    bool coupled = td.cloud().solution().coupled();
-
     // check if parcel belongs to liquid core
     if (liquidCore() > 0.5)
     {
-        // liquid core parcels should not interact with the gas
-        td.cloud().solution().coupled() = false;
+        // liquid core parcels should not experience coupled forces
+        td.cloud().forces().setCalcCoupled(false);
     }
 
     // get old mixture composition
@@ -138,8 +136,8 @@ void Foam::SprayParcel<ParcelType>::calc
         }
     }
 
-    // restore coupled
-    td.cloud().solution().coupled() = coupled;
+    // restore coupled forces
+    td.cloud().forces().setCalcCoupled(true);
 }
 
 
@@ -246,7 +244,7 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
     const scalar mass = p.mass();
     const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
     const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
-    scalar tMom = mass/(Fcp.Sp() + Fncp.Sp());
+    this->tMom() = mass/(Fcp.Sp() + Fncp.Sp());
 
     const vector g = td.cloud().g().value();
 
@@ -274,7 +272,7 @@ void Foam::SprayParcel<ParcelType>::calcBreakup
             muAv,
             Urel,
             Urmag,
-            tMom,
+            this->tMom(),
             dChild,
             massChild
         )
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C b/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C
index c342455136..63b5d4ff5c 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C
+++ b/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C
@@ -48,7 +48,7 @@ Foam::BlobsSheetAtomization<CloudType>::BlobsSheetAtomization
 :
     AtomizationModel<CloudType>(am),
     B_(am.B_),
-    angle_(am.B_)
+    angle_(am.angle_)
 {}
 
 
diff --git a/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H b/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H
index 0927612564..22970e7c48 100644
--- a/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H
+++ b/src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -112,7 +112,7 @@ public:
             const scalar rho,
             const scalar mu,
             const scalar sigma,
-            const scalar volHlowRate,
+            const scalar volFlowRate,
             const scalar rhoAv,
             const scalar Urel,
             const vector& pos,
diff --git a/src/meshTools/coordinateSystems/coordinateRotation/localAxesRotation.C b/src/meshTools/coordinateSystems/coordinateRotation/localAxesRotation.C
index 305601f116..719026227e 100644
--- a/src/meshTools/coordinateSystems/coordinateRotation/localAxesRotation.C
+++ b/src/meshTools/coordinateSystems/coordinateRotation/localAxesRotation.C
@@ -52,7 +52,8 @@ namespace Foam
 
 Foam::localAxesRotation::localAxesRotation
 (
-    const dictionary& dict, const objectRegistry& orb
+    const dictionary& dict,
+    const objectRegistry& orb
 )
 :
     Rptr_(),
@@ -88,10 +89,10 @@ Foam::localAxesRotation::localAxesRotation
     e3_()
 {
     FatalErrorIn("localAxesRotation(const dictionary&)")
-        << " localAxesRotation can not be contructed from  dictionary "
+        << " localAxesRotation can not be constructed from  dictionary "
         << " use the construtctor : "
            "("
-           "    const dictionary& dict, const objectRegistry& orb"
+           "    const dictionary&, const objectRegistry&"
            ")"
         << exit(FatalIOError);
 }
@@ -112,7 +113,7 @@ Foam::vector Foam::localAxesRotation::transform(const vector& st) const
 {
     notImplemented
     (
-        "vector Foam::localAxesRotation::transform(const vector&) const"
+        "vector localAxesRotation::transform(const vector&) const"
     );
     return vector::zero;
 }
@@ -122,7 +123,7 @@ Foam::vector Foam::localAxesRotation::invTransform(const vector& st) const
 {
     notImplemented
     (
-        "vector Foam::localAxesRotation::invTransform(const vector&) const"
+        "vector localAxesRotation::invTransform(const vector&) const"
     );
     return vector::zero;
 }
@@ -135,7 +136,10 @@ Foam::tmp<Foam::vectorField> Foam::localAxesRotation::transform
 {
     if (Rptr_->size() != st.size())
     {
-        FatalErrorIn("localAxesRotation::transform(const vectorField&)")
+        FatalErrorIn
+        (
+            "tmp<vectorField> localAxesRotation::transform(const vectorField&)"
+        )
             << "vectorField st has different size to tensorField "
             << abort(FatalError);
     }
@@ -160,7 +164,13 @@ Foam::tmp<Foam::tensorField> Foam::localAxesRotation::transformTensor
 {
     if (Rptr_->size() != st.size())
     {
-        FatalErrorIn("localAxesRotation::transformTensor(const tensorField&)")
+        FatalErrorIn
+        (
+            "tmp<tensorField> localAxesRotation::transformTensor"
+            "("
+                "const tensorField&"
+            ")"
+        )
             << "tensorField st has different size to tensorField Tr"
             << abort(FatalError);
     }
@@ -173,7 +183,11 @@ Foam::tensor Foam::localAxesRotation::transformTensor
     const tensor& st
 ) const
 {
-    notImplemented("tensor localAxesRotation::transformTensor() const");
+    notImplemented
+    (
+        "tensor localAxesRotation::transformTensor(const tensor&) const"
+    );
+
     return tensor::zero;
 }
 
@@ -186,7 +200,14 @@ Foam::tmp<Foam::tensorField> Foam::localAxesRotation::transformTensor
 {
     if (cellMap.size() != st.size())
     {
-        FatalErrorIn("localAxesRotation::transformTensor(const tensorField&)")
+        FatalErrorIn
+        (
+            "tmp<tensorField> localAxesRotation::transformTensor"
+            "("
+                "const tensorField&"
+                "const labelList&"
+            ")"
+        )
             << "tensorField st has different size to tensorField Tr"
             << abort(FatalError);
     }
diff --git a/src/meshTools/searchableSurface/triSurfaceMesh.H b/src/meshTools/searchableSurface/triSurfaceMesh.H
index 150bd768c1..5751ead958 100644
--- a/src/meshTools/searchableSurface/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurface/triSurfaceMesh.H
@@ -85,6 +85,7 @@ private:
         //- Is surface closed
         mutable label surfaceClosed_;
 
+
     // Private Member Functions
 
         ////- Helper: find instance of files without header
@@ -166,7 +167,7 @@ public:
         //- Move points
         virtual void movePoints(const pointField&);
 
-        //- Demand driven contruction of octree for boundary edges
+        //- Demand driven construction of octree for boundary edges
         const indexedOctree<treeDataEdge>& edgeTree() const;
 
 
@@ -260,6 +261,7 @@ public:
                 List<volumeType>&
             ) const;
 
+
         // Other
 
             //- WIP. Store element-wise field.
@@ -285,7 +287,6 @@ public:
                 IOstream::versionNumber ver,
                 IOstream::compressionType cmp
             ) const;
-
 };
 
 
diff --git a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/patchifyObstacles b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/patchifyObstacles
index bc9fe0a8a3..f9218e997d 100755
--- a/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/patchifyObstacles
+++ b/tutorials/lagrangian/reactingParcelFilmFoam/hotBoxes/patchifyObstacles
@@ -155,7 +155,7 @@ createSetsAndZone floor -$tol $floorMax -$tol $floorMax -$tol $tol
 echo "cellSet floorCells new faceToCell floorFaces owner" >> $tmpSetSet
 echo "faceZoneSet floorFaces new setsToFaceZone floorFaces floorCells" >> $tmpSetSet
 
-setSet -batch $tmpSetSet > log.setSet.patchifyObstacles >/dev/null 2>&1
+setSet -batch $tmpSetSet > log.setSet.patchifyObstacles 2>&1
 
 
 # *************************************************************************