diff --git a/applications/utilities/mesh/generation/foamyHexMesh/Make/options b/applications/utilities/mesh/generation/foamyHexMesh/Make/options
index 9461fa3725d8b76cc7cdce5cc102b5aa1933bdbf..f03dbf90976d5f7ecedf098b24beced6e5cb968c 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/Make/options
+++ b/applications/utilities/mesh/generation/foamyHexMesh/Make/options
@@ -21,6 +21,7 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
     -IvectorTools
 
 EXE_LIBS = \
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C
index 551380c6e2017e2c7256482d8255ff63e1980f76..1b570e90e768b204555113f8e436687bad288563 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMeshTools/DelaunayMeshToolsTemplates.C
@@ -279,11 +279,9 @@ Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
     const Triangulation& t
 )
 {
-    tmp<pointField> tpts(new pointField(t.number_of_vertices(), point::max));
+    tmp<pointField> tpts(new pointField(t.vertexCount(), point::max));
     pointField& pts = tpts();
 
-    label nVert = 0;
-
     for
     (
         typename Triangulation::Finite_vertices_iterator vit =
@@ -292,9 +290,9 @@ Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
         ++vit
     )
     {
-        if (vit->internalOrBoundaryPoint())
+        if (vit->internalOrBoundaryPoint() && !vit->referred())
         {
-            pts[nVert++] = topoint(vit->point());
+            pts[vit->index()] = topoint(vit->point());
         }
     }
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options
index 9963458aa5b3518a505eeb095fe0874c7105c3fa..2be5855323b10652c1de0a434a8d1ec95f6f2c16 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options
@@ -22,6 +22,7 @@ EXE_INC = \
     -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
     -IPrintTable \
     -I../vectorTools
 
@@ -32,4 +33,5 @@ LIB_LIBS = \
     -ltriSurface \
     -ldynamicMesh \
     -lsurfMesh \
-    -lsampling
+    -lsampling \
+    -lautoMesh
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index dcaf44f99754f758b63723aa61180e140fe3905e..6dd0675993f1835726ab08988eb01e56afd949f0 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -600,6 +600,29 @@ private:
             PackedBoolList& boundaryFacesToRemove
         );
 
+        //- From meshRefinementBaffles.C. Use insidePoint for a surface to
+        //  determine the cell zone.
+        void findCellZoneInsideWalk
+        (
+            const polyMesh& mesh,
+            const labelList& locationSurfaces,
+            const labelList& faceToSurface,
+            labelList& cellToSurface
+        ) const;
+
+        //- Calculate the cell zones from cellCentres using all closed surfaces
+        labelList calcCellZones(const pointField& cellCentres) const;
+
+        //- Calculate the face zones
+        void calcFaceZones
+        (
+            const polyMesh& mesh,
+            const pointField& cellCentres,
+            const labelList& cellToSurface,
+            labelList& faceToSurface,
+            boolList& flipMap
+        ) const;
+
         //- Tet mesh calculation
         void calcTetMesh
         (
@@ -712,9 +735,6 @@ private:
             bool includeEmptyPatches = false
         ) const;
 
-        //- Create the cell centres to use for the mesh
-        void createCellCentres(pointField& cellCentres) const;
-
         //- Sort the faces, owner and neighbour lists into
         //  upper-triangular order.  For internal faces only, use
         //  before adding patch faces
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index 86892c28d940c822f80756da75be935159891ce8..c943ca978b9e5abb538a8febb19a2040bfafc3a8 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -2852,32 +2852,6 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 }
 
 
-void Foam::conformalVoronoiMesh::createCellCentres
-(
-    pointField& cellCentres
-) const
-{
-    cellCentres.setSize(number_of_vertices(), point::max);
-
-    label vertI = 0;
-
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (vit->internalOrBoundaryPoint())
-        {
-            cellCentres[vertI++] = topoint(vit->point());
-        }
-    }
-
-    cellCentres.setSize(vertI);
-}
-
-
 void Foam::conformalVoronoiMesh::sortFaces
 (
     faceList& faces,
@@ -3095,7 +3069,7 @@ Foam::labelList Foam::conformalVoronoiMesh::removeUnusedCells
 {
     Info<< nl << "Removing unused cells" << endl;
 
-    PackedBoolList cellUsed(number_of_vertices(), false);
+    PackedBoolList cellUsed(vertexCount(), false);
 
     // Scan all faces to find all of the cells that are used
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index bb34cbd4a69c0e92987c75c77af047b55dd009a5..07fd509e1f5e100151f2eaa9fff5cda5b26e0585 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -35,6 +35,11 @@ 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"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -402,6 +407,376 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
 }
 
 
+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(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();
+
+    forAll(faces, faceI)
+    {
+        const label ownerSurfaceI = cellToSurface[faceOwner[faceI]];
+
+        if (mesh.isInternalFace(faceI))
+        {
+            const label neiSurfaceI = cellToSurface[faceNeighbour[faceI]];
+
+            flipMap[faceI] =
+                (
+                    ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
+                  ? false
+                  : true
+                );
+
+            if
+            (
+                (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
+             && ownerSurfaceI != neiSurfaceI
+            )
+            {
+                if (ownerSurfaceI > neiSurfaceI)
+                {
+                    faceToSurface[faceI] = ownerSurfaceI;
+                }
+                else
+                {
+                    faceToSurface[faceI] = neiSurfaceI;
+                }
+            }
+        }
+        else
+        {
+            if (ownerSurfaceI >= 0)
+            {
+                faceToSurface[faceI] = ownerSurfaceI;
+            }
+        }
+    }
+
+
+    const PtrList<surfaceZonesInfo>& surfZones =
+        geometryToConformTo().surfZones();
+
+    labelList insidePointNamedSurfaces
+    (
+        surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
+    );
+
+    // Use intersection of cellCentre connections
+    forAll(faces, faceI)
+    {
+        if
+        (
+            mesh.isInternalFace(faceI)
+         && faceToSurface[faceI] < 0
+        )
+        {
+            const label own = faceOwner[faceI];
+            const label nei = faceNeighbour[faceI];
+
+            pointIndexHit surfHit;
+            label hitSurface;
+
+            geometryToConformTo().findSurfaceAnyIntersection
+            (
+                cellCentres[own],
+                cellCentres[nei],
+                surfHit,
+                hitSurface
+            );
+
+            if (surfHit.hit())
+            {
+                if (findIndex(insidePointNamedSurfaces, hitSurface) != -1)
+                {
+                    faceToSurface[faceI] = hitSurface;
+
+                    vectorField norm;
+                    geometryToConformTo().getNormal
+                    (
+                        hitSurface,
+                        List<pointIndexHit>(1, surfHit),
+                        norm
+                    );
+
+                    const vector fN = faces[faceI].normal(mesh.points());
+
+                    if ((norm[0] & fN/(mag(fN) + SMALL)) < 0)
+                    {
+                        flipMap[faceI] = true;
+                    }
+                    else
+                    {
+                        flipMap[faceI] = false;
+                    }
+                }
+            }
+        }
+    }
+
+
+    labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    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,
@@ -883,7 +1258,11 @@ void Foam::conformalVoronoiMesh::writeMesh
             // Check that the patch is not empty on every processor
             reduce(totalPatchSize, sumOp<label>());
 
-            if (totalPatchSize > 0)
+            if
+            (
+                totalPatchSize > 0
+//             && !geometryToConformTo().surfZones().set(p)
+            )
             {
                 patches[nValidPatches] = polyPatch::New
                 (
@@ -898,6 +1277,145 @@ void Foam::conformalVoronoiMesh::writeMesh
         }
     }
 
+    patches.setSize(nValidPatches);
+
+    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)
+        {
+            if (!mesh.isInternalFace(faceI))
+            {
+                continue;
+            }
+
+            label surfaceI = faceToSurface[faceI];
+
+            if (surfaceI >= 0)
+            {
+                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
+                    )
+                );
+            }
+        }
+
+        // Change the mesh (no inflation, parallel sync)
+        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
+    }
+
+
+
+
     // Add indirectPatchFaces to a face zone
     {
         labelList addr(boundaryFacesToRemove.count());
@@ -928,10 +1446,6 @@ void Foam::conformalVoronoiMesh::writeMesh
         );
     }
 
-    patches.setSize(nValidPatches);
-
-    mesh.addFvPatches(patches);
-
     timeCheck("Before fvMesh filtering");
 
     autoPtr<polyMeshFilter> meshFilter;
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 0e4b857d4dab316a880baeeca6ec1fe581a1155f..351b85c952dd13a9f9875c3e6a86f3c346ea8c57 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -579,30 +579,6 @@ Foam::conformationSurfaces::conformationSurfaces
             Info<< features_[fI].name() << endl;
         }
     }
-
-    Info<< "ZONES" << endl;
-    forAll(surfZones_, surfI)
-    {
-        if (surfZones_.set(surfI))
-        {
-            const surfaceZonesInfo& sInfo = surfZones_[surfI];
-
-            Info<< "    " << surfI << nl
-                << "    faceZoneName    = " << sInfo.faceZoneName() << nl
-                << "    cellZoneName    = " << sInfo.cellZoneName() << nl
-                << "    zoneInside      = "
-                << surfaceZonesInfo::areaSelectionAlgoNames[sInfo.zoneInside()]
-                << nl
-                << "    zoneInsidePoint = " << sInfo.zoneInsidePoint() << nl
-                << "    faceType        = "
-                << surfaceZonesInfo::faceZoneTypeNames[sInfo.faceType()]
-                << endl;
-        }
-        else
-        {
-            Info<< "    " << surfI << " EMPTY" << endl;
-        }
-    }
 }
 
 
diff --git a/applications/utilities/mesh/generation/foamyQuadMesh/Make/options b/applications/utilities/mesh/generation/foamyQuadMesh/Make/options
index d1a5bfde5422911eb54e6feed45f9fbcd99a6862..62aec78cdb2e485624764ac37477d59dd9387207 100755
--- a/applications/utilities/mesh/generation/foamyQuadMesh/Make/options
+++ b/applications/utilities/mesh/generation/foamyQuadMesh/Make/options
@@ -23,6 +23,7 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/mesh/autoMesh/lnInclude
 
 EXE_LIBS = \
     $(CGAL_LIBS) \
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
index 4a8b82e129b0da9215d41fd7d56cdb9c8bc00e68..5d2b1c2cfcceb84fc560bc238c6496562ebe5350 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
@@ -124,7 +124,7 @@ Foam::surfaceZonesInfo::surfaceZonesInfo
             {
                 IOWarningIn
                 (
-                    "refinementSurfaces::refinementSurfaces(..)",
+                    "surfaceZonesInfo::surfaceZonesInfo(..)",
                     surfacesDict
                 )   << "Illegal entry zoneInside "
                     << areaSelectionAlgoNames[zoneInside_]
@@ -137,7 +137,7 @@ Foam::surfaceZonesInfo::surfaceZonesInfo
         {
             IOWarningIn
             (
-                "refinementSurfaces::refinementSurfaces(..)",
+                "surfaceZonesInfo::surfaceZonesInfo(..)",
                 surfacesDict
             )   << "Unused entry zoneInside for faceZone "
                 << faceZoneName_
@@ -215,7 +215,11 @@ Foam::labelList Foam::surfaceZonesInfo::getNamedSurfaces
     label namedI = 0;
     forAll(surfList, surfI)
     {
-        if (surfList[surfI].faceZoneName().size())
+        if
+        (
+            surfList.set(surfI)
+         && surfList[surfI].faceZoneName().size()
+        )
         {
             namedSurfaces[namedI++] = surfI;
         }
@@ -241,7 +245,8 @@ Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
     {
         if
         (
-            surfList[surfI].cellZoneName().size()
+            surfList.set(surfI)
+         && surfList[surfI].cellZoneName().size()
          && (
                 surfList[surfI].zoneInside() == surfaceZonesInfo::INSIDE
              || surfList[surfI].zoneInside() == surfaceZonesInfo::OUTSIDE
@@ -258,6 +263,35 @@ Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
 }
 
 
+// Get indices of closed named surfaces
+Foam::labelList Foam::surfaceZonesInfo::getAllClosedNamedSurfaces
+(
+    const PtrList<surfaceZonesInfo>& surfList,
+    const searchableSurfaces& allGeometry,
+    const labelList& surfaces
+)
+{
+    labelList closed(surfList.size());
+
+    label closedI = 0;
+    forAll(surfList, surfI)
+    {
+        if
+        (
+            surfList.set(surfI)
+         && surfList[surfI].cellZoneName().size()
+         && allGeometry[surfaces[surfI]].hasVolumeType()
+        )
+        {
+            closed[closedI++] = surfI;
+        }
+    }
+    closed.setSize(closedI);
+
+    return closed;
+}
+
+
 // Get indices of named surfaces with a
 Foam::labelList Foam::surfaceZonesInfo::getInsidePointNamedSurfaces
 (
@@ -271,7 +305,8 @@ Foam::labelList Foam::surfaceZonesInfo::getInsidePointNamedSurfaces
     {
         if
         (
-            surfList[surfI].cellZoneName().size()
+            surfList.set(surfI)
+         && surfList[surfI].cellZoneName().size()
          && surfList[surfI].zoneInside() == surfaceZonesInfo::INSIDEPOINT
         )
         {
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H
index dbfd86e9f855d79fcb3a6816b6e1394147370f4b..39a48adcfd8c33983c21633ac439a1cbc2265cda 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H
@@ -199,6 +199,14 @@ public:
                 const labelList& surfaces
             );
 
+            //- Get indices of surfaces with a cellZone that are closed.
+            static labelList getAllClosedNamedSurfaces
+            (
+                const PtrList<surfaceZonesInfo>& surfList,
+                const searchableSurfaces& allGeometry,
+                const labelList& surfaces
+            );
+
             //- Get indices of surfaces with a cellZone that have 'insidePoint'
             //  section.
             static labelList getInsidePointNamedSurfaces