diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 6dd0675993f1835726ab08988eb01e56afd949f0..9a22055dc90d708599b5a0573f6a2d59fdbb37aa 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -600,6 +600,19 @@ private:
             PackedBoolList& boundaryFacesToRemove
         );
 
+        void calcNeighbourCellCentres
+        (
+            const polyMesh& mesh,
+            const pointField& cellCentres,
+            pointField& neiCc
+        ) const;
+
+        void selectSeparatedCoupledFaces
+        (
+            const polyMesh& mesh,
+            boolList& selected
+        ) const;
+
         //- From meshRefinementBaffles.C. Use insidePoint for a surface to
         //  determine the cell zone.
         void findCellZoneInsideWalk
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index c943ca978b9e5abb538a8febb19a2040bfafc3a8..4e70aab9649b6ea57fc5c06dc494f1008e5d6cf9 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -2744,12 +2744,45 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                 }
                 else
                 {
-                    // internal face
-                    faces[dualFaceI] = newDualFace;
-                    owner[dualFaceI] = own;
-                    neighbour[dualFaceI] = nei;
+//                    if
+//                    (
+//                        ptPairs_.isPointPair(vA, vB)
+//                     || ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
+//                    )
+//                    {
+                        patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
+//                    }
+
+                    if (patchIndex != -1)
+                    {
+//                        patchFaces[patchIndex].append(newDualFace);
+//                        patchOwners[patchIndex].append(own);
+//                        indirectPatchFace[patchIndex].append(false);
+//
+//                        if
+//                        (
+//                            labelPair(vB->index(), vB->procIndex())
+//                          < labelPair(vA->index(), vA->procIndex())
+//                        )
+//                        {
+//                            patchPPSlaves[patchIndex].append(vB->index());
+//                        }
+//                        else
+//                        {
+//                            patchPPSlaves[patchIndex].append(vA->index());
+//                        }
 
-                    dualFaceI++;
+//                        baffleFaces[dualFaceI] = patchIndex;
+                    }
+//                    else
+                    {
+                        // internal face
+                        faces[dualFaceI] = newDualFace;
+                        owner[dualFaceI] = own;
+                        neighbour[dualFaceI] = nei;
+
+                        dualFaceI++;
+                    }
                 }
             }
         }
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 07fd509e1f5e100151f2eaa9fff5cda5b26e0585..dc272a89d6f2eb8f6c4077b36d9f1ea52e164367 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -40,6 +40,7 @@ License
 #include "polyModifyFace.H"
 #include "syncTools.H"
 #include "regionSplit.H"
+#include "OBJstream.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -111,7 +112,7 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
 
     {
         pointField points;
-        labelList boundaryPts(number_of_finite_cells(), -1);
+        labelList boundaryPts;
         faceList faces;
         labelList owner;
         labelList neighbour;
@@ -407,6 +408,78 @@ 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,
@@ -417,7 +490,7 @@ void Foam::conformalVoronoiMesh::findCellZoneInsideWalk
 {
     // Analyse regions. Reuse regionsplit
     boolList blockedFace(mesh.nFaces());
-    //selectSeparatedCoupledFaces(blockedFace);
+    selectSeparatedCoupledFaces(mesh, blockedFace);
 
     forAll(faceToSurface, faceI)
     {
@@ -630,42 +703,90 @@ void Foam::conformalVoronoiMesh::calcFaceZones
     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]];
 
-            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;
-                }
+                flipMap[faceI] =
+                    (
+                        ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
+                      ? false
+                      : true
+                    );
+
+                faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
             }
         }
         else
         {
-            if (ownerSurfaceI >= 0)
+            label patchID = mesh.boundaryMesh().whichPatch(faceI);
+
+            if (mesh.boundaryMesh()[patchID].coupled())
             {
-                faceToSurface[faceI] = ownerSurfaceI;
+                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;
+                }
             }
         }
     }
@@ -674,104 +795,151 @@ void Foam::conformalVoronoiMesh::calcFaceZones
     const PtrList<surfaceZonesInfo>& surfZones =
         geometryToConformTo().surfZones();
 
-    labelList insidePointNamedSurfaces
+    labelList unclosedSurfaces
     (
-        surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
+        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
-        (
-            mesh.isInternalFace(faceI)
-         && faceToSurface[faceI] < 0
-        )
+        if (faceToSurface[faceI] >= 0)
         {
-            const label own = faceOwner[faceI];
-            const label nei = faceNeighbour[faceI];
+            continue;
+        }
 
-            pointIndexHit surfHit;
-            label hitSurface;
+        label patchID = mesh.boundaryMesh().whichPatch(faceI);
 
-            geometryToConformTo().findSurfaceAnyIntersection
+        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.hit())
+            if (surfHit.size() == 1 && surfHit[0].hit())
             {
-                if (findIndex(insidePointNamedSurfaces, hitSurface) != -1)
-                {
-                    faceToSurface[faceI] = hitSurface;
-
-                    vectorField norm;
-                    geometryToConformTo().getNormal
+                intersections.write
+                (
+                    linePointRef
                     (
-                        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;
-                    }
-                }
+                        cellCentres[own],
+                        neiCc[faceI - mesh.nInternalFaces()]
+                    )
+                );
             }
         }
-    }
-
 
-    labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
+        // If there are multiple intersections then do not add to
+        // a faceZone
+        if (surfHit.size() == 1 && surfHit[0].hit())
         {
-            forAll(pp, i)
+            if (findIndex(unclosedSurfaces, hitSurface[0]) != -1)
             {
-                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];
+                vectorField norm;
+                geometryToConformTo().getNormal
+                (
+                    hitSurface[0],
+                    List<pointIndexHit>(1, surfHit[0]),
+                    norm
+                );
 
-        if (pp.coupled())
-        {
-            forAll(pp, i)
-            {
-                label faceI = pp.start()+i;
-                label ownSurface = cellToSurface[faceOwner[faceI]];
-                label neiSurface = neiCellSurface[faceI-mesh.nInternalFaces()];
+                vector fN = faces[faceI].normal(mesh.points());
+                fN /= mag(fN) + SMALL;
 
-                if (faceToSurface[faceI] == -1 && (ownSurface != neiSurface))
+                if ((norm[0] & fN) < 0)
                 {
-                    // Give face the max cell zone
-                    faceToSurface[faceI] =  max(ownSurface, neiSurface);
+                    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>());
 }
@@ -1049,7 +1217,7 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
 
     pBufs.finishedSends();
 
-    Info<< incrIndent << indent << "    Face ordering initialised..." << endl;
+    Info<< incrIndent << indent << "Face ordering initialised..." << endl;
 
     // Receive and calculate ordering
     bool anyChanged = false;
@@ -1108,7 +1276,7 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
         }
     }
 
-    Info<< incrIndent << indent << "    Faces matched." << endl;
+    Info<< incrIndent << indent << "Faces matched." << endl;
 
     reduce(anyChanged, orOp<bool>());
 
@@ -1145,6 +1313,7 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
             << " faces have been reordered" << nl
             << indent << returnReduce(nRotated, sumOp<label>())
             << " faces have been rotated"
+            << decrIndent << decrIndent
             << decrIndent << decrIndent << endl;
     }
 }
@@ -1192,7 +1361,8 @@ void Foam::conformalVoronoiMesh::writeMesh
         );
     }
 
-    Info<< "    Constructing mesh" << endl;
+    Info<< incrIndent;
+    Info<< indent << "Constructing mesh" << endl;
 
     timeCheck("Before fvMesh construction");
 
@@ -1212,7 +1382,7 @@ void Foam::conformalVoronoiMesh::writeMesh
         xferMove(neighbour)
     );
 
-    Info<< "    Adding patches to mesh" << endl;
+    Info<< indent << "Adding patches to mesh" << endl;
 
     List<polyPatch*> patches(patchNames.size());
 
@@ -1258,11 +1428,7 @@ void Foam::conformalVoronoiMesh::writeMesh
             // Check that the patch is not empty on every processor
             reduce(totalPatchSize, sumOp<label>());
 
-            if
-            (
-                totalPatchSize > 0
-//             && !geometryToConformTo().surfZones().set(p)
-            )
+            if (totalPatchSize > 0)
             {
                 patches[nValidPatches] = polyPatch::New
                 (
@@ -1379,14 +1545,16 @@ void Foam::conformalVoronoiMesh::writeMesh
 
         forAll(faceToSurface, faceI)
         {
-            if (!mesh.isInternalFace(faceI))
+            label surfaceI = faceToSurface[faceI];
+
+            if (surfaceI < 0)
             {
                 continue;
             }
 
-            label surfaceI = faceToSurface[faceI];
+            label patchID = mesh.boundaryMesh().whichPatch(faceI);
 
-            if (surfaceI >= 0)
+            if (mesh.isInternalFace(faceI))
             {
                 label own = faceOwner[faceI];
                 label nei = faceNeighbour[faceI];
@@ -1407,6 +1575,26 @@ void Foam::conformalVoronoiMesh::writeMesh
                     )
                 );
             }
+            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)
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 351b85c952dd13a9f9875c3e6a86f3c346ea8c57..43d9c770134a330f62629b67c488b32a91a1173d 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -872,6 +872,42 @@ void Foam::conformationSurfaces::findSurfaceAnyIntersection
 }
 
 
+void Foam::conformationSurfaces::findSurfaceAllIntersections
+(
+    const point& start,
+    const point& end,
+    List<pointIndexHit>& surfHit,
+    labelList& hitSurface
+) const
+{
+    labelListList hitSurfaces;
+    List<List<pointIndexHit> > hitInfo;
+
+    searchableSurfacesQueries::findAllIntersections
+    (
+        allGeometry_,
+        surfaces_,
+        pointField(1, start),
+        pointField(1, end),
+        hitSurfaces,
+        hitInfo
+    );
+
+    surfHit = hitInfo[0];
+
+    hitSurface.setSize(hitSurfaces[0].size());
+
+    forAll(hitSurfaces[0], surfI)
+    {
+        // hitSurfaces has returned the index of the entry in surfaces_ that was
+        // found, not the index of the surface in allGeometry_, translating this
+        // to allGeometry_
+
+        hitSurface[surfI] = surfaces_[hitSurfaces[0][surfI]];
+    }
+}
+
+
 void Foam::conformationSurfaces::findSurfaceNearestIntersection
 (
     const point& start,
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index 68c5fd28b92e93c159622fa21ca2883da222cbd0..12932255da25a69e5effc0b1bbd81c223dac35be 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -257,6 +257,14 @@ public:
                 label& hitSurface
             ) const;
 
+            void findSurfaceAllIntersections
+            (
+                const point& start,
+                const point& end,
+                List<pointIndexHit>& surfHit,
+                labelList& hitSurface
+            ) const;
+
             //- Finding the nearestIntersection of the surface to start
             void findSurfaceNearestIntersection
             (
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index b52777ffb468ba44de64d0422e208db2ddab3d9c..a83b2bf1909f9e962ac843bd95e981d290148a92 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -847,7 +847,9 @@ Foam::label Foam::globalMeshData::findTransform
     {
         FatalErrorIn("globalMeshData::findTransform(..)")
             << "Problem. Cannot find " << remotePoint
-            << " or " << localPoint << " in " << info
+            << " or " << localPoint  << " "
+            << coupledPatch().localPoints()[localPoint]
+            << " in " << info
             << endl
             << "remoteTransformI:" << remoteTransformI << endl
             << "localTransformI:" << localTransformI
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
index 5d2b1c2cfcceb84fc560bc238c6496562ebe5350..69c3f7777b95ea4b1a02a1d3d817673d55a7779b 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
@@ -182,7 +182,6 @@ Foam::surfaceZonesInfo::surfaceZonesInfo(const surfaceZonesInfo& surfZone)
 {}
 
 
-// Get indices of unnamed surfaces (surfaces without faceZoneName)
 Foam::labelList Foam::surfaceZonesInfo::getUnnamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList
@@ -204,7 +203,6 @@ Foam::labelList Foam::surfaceZonesInfo::getUnnamedSurfaces
 }
 
 
-// Get indices of named surfaces (surfaces with faceZoneName)
 Foam::labelList Foam::surfaceZonesInfo::getNamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList
@@ -230,7 +228,6 @@ Foam::labelList Foam::surfaceZonesInfo::getNamedSurfaces
 }
 
 
-// Get indices of closed named surfaces
 Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList,
@@ -263,7 +260,33 @@ Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
 }
 
 
-// Get indices of closed named surfaces
+Foam::labelList Foam::surfaceZonesInfo::getUnclosedNamedSurfaces
+(
+    const PtrList<surfaceZonesInfo>& surfList,
+    const searchableSurfaces& allGeometry,
+    const labelList& surfaces
+)
+{
+    labelList unclosed(surfList.size());
+
+    label unclosedI = 0;
+    forAll(surfList, surfI)
+    {
+        if
+        (
+            surfList.set(surfI)
+         && !allGeometry[surfaces[surfI]].hasVolumeType()
+        )
+        {
+            unclosed[unclosedI++] = surfI;
+        }
+    }
+    unclosed.setSize(unclosedI);
+
+    return unclosed;
+}
+
+
 Foam::labelList Foam::surfaceZonesInfo::getAllClosedNamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList,
@@ -292,7 +315,6 @@ Foam::labelList Foam::surfaceZonesInfo::getAllClosedNamedSurfaces
 }
 
 
-// Get indices of named surfaces with a
 Foam::labelList Foam::surfaceZonesInfo::getInsidePointNamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H
index 39a48adcfd8c33983c21633ac439a1cbc2265cda..eb59aad34fab72539af60b17a02d55a2c05051e2 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 unclosed
+            static labelList getUnclosedNamedSurfaces
+            (
+                const PtrList<surfaceZonesInfo>& surfList,
+                const searchableSurfaces& allGeometry,
+                const labelList& surfaces
+            );
+
             //- Get indices of surfaces with a cellZone that are closed.
             static labelList getAllClosedNamedSurfaces
             (