From 615a74626e15e63daf8be2769a0a2c9b65943160 Mon Sep 17 00:00:00 2001
From: laurence <laurence>
Date: Wed, 25 Sep 2013 12:37:19 +0100
Subject: [PATCH] ENH: foamyHexMesh: Include baffle handling

---
 .../advanced/collapseEdges/collapseEdges.C    |  45 +-
 .../mesh/generation/foamyHexMesh/Make/options |   1 +
 .../DelaunayMesh/DelaunayMesh.C               |  29 +-
 .../DelaunayMesh/DelaunayMesh.H               |   8 +-
 .../DelaunayMesh/DistributedDelaunayMesh.C    |   3 +-
 .../conformalVoronoiMesh/Make/files           |   1 +
 .../conformalVoronoiMesh/Make/options         |   1 +
 .../cellShapeControlMesh.C                    | 117 ++--
 .../cellShapeControlMesh.H                    |  11 +-
 .../searchableSurfaceControl.C                |   5 +
 .../controlMeshRefinement.C                   |   2 +-
 .../conformalVoronoiMesh.C                    | 173 +++--
 .../conformalVoronoiMesh.H                    |  48 +-
 .../conformalVoronoiMeshCalcDualMesh.C        | 607 ++----------------
 .../conformalVoronoiMeshConformToSurface.C    | 452 ++++++-------
 .../conformalVoronoiMeshFeaturePoints.C       | 443 +++++++++++--
 .../conformalVoronoiMeshI.H                   | 109 +---
 .../conformalVoronoiMeshIO.C                  | 142 +++-
 .../featurePointConformer.C                   |  12 +-
 .../featurePointConformer.H                   |   7 +
 .../featurePointConformerI.H                  |   6 +
 .../featurePointConformerSpecialisations.C    | 196 +++++-
 .../indexedCell/indexedCell.H                 |   2 +
 .../indexedCell/indexedCellI.H                |  17 +
 .../conformationSurfaces.C                    |  29 +-
 .../autoDensity/autoDensity.C                 |  20 +-
 .../pointPairs/pointPairs.C                   | 244 +++++++
 .../pointPairs/pointPairs.H                   | 163 +++++
 .../searchableBoxFeatures.C                   |  11 +-
 .../searchablePlateFeatures.C                 | 218 +++++++
 .../searchablePlateFeatures.H                 | 120 ++++
 .../triSurfaceMeshFeatures.C                  |  24 +-
 .../triSurfaceMeshFeatures.H                  |   3 +
 .../surfaceBooleanFeatures.C                  | 285 ++++++--
 .../surfaceFeatureExtract.C                   |  29 +-
 .../extendedFeatureEdgeMesh.C                 |  32 +-
 .../extendedFeatureEdgeMesh.H                 |  31 +-
 .../extendedFeatureEdgeMeshI.H                |  23 +
 .../surfaceIntersection/edgeIntersections.C   | 136 ++--
 .../surfaceIntersection/surfaceIntersection.C |   7 +-
 .../surfaceIntersection/surfaceIntersection.H |   2 +-
 .../mixerVessel/system/foamyHexMeshDict       |  10 +-
 42 files changed, 2444 insertions(+), 1380 deletions(-)
 create mode 100644 applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.C
 create mode 100644 applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.H
 create mode 100644 applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.C
 create mode 100644 applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.H

diff --git a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
index 0515205c78..5b24606613 100644
--- a/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
+++ b/applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C
@@ -52,6 +52,7 @@ Usage
 #include "polyTopoChange.H"
 #include "fvMesh.H"
 #include "polyMeshFilter.H"
+#include "faceSet.H"
 
 using namespace Foam;
 
@@ -74,9 +75,9 @@ int main(int argc, char *argv[])
 
     argList::addOption
     (
-        "collapseFaceZone",
-        "zoneName",
-        "Collapse faces that are in the supplied face zone"
+        "collapseFaceSet",
+        "faceSet",
+        "Collapse faces that are in the supplied face set"
     );
 
 #   include "addOverwriteOption.H"
@@ -93,15 +94,15 @@ int main(int argc, char *argv[])
     const bool overwrite = args.optionFound("overwrite");
 
     const bool collapseFaces = args.optionFound("collapseFaces");
-    const bool collapseFaceZone = args.optionFound("collapseFaceZone");
+    const bool collapseFaceSet = args.optionFound("collapseFaceSet");
 
-    if (collapseFaces && collapseFaceZone)
+    if (collapseFaces && collapseFaceSet)
     {
         FatalErrorIn("main(int, char*[])")
             << "Both face zone collapsing and face collapsing have been"
             << "selected. Choose only one of:" << nl
             << "    -collapseFaces" << nl
-            << "    -collapseFaceZone <faceZoneName>"
+            << "    -collapseFaceSet <faceSet>"
             << abort(FatalError);
     }
 
@@ -117,7 +118,6 @@ int main(int argc, char *argv[])
         ),
         labelList(mesh.nPoints(), labelMin)
     );
-
     forAll(timeDirs, timeI)
     {
         runTime.setTime(timeDirs[timeI], timeI);
@@ -128,6 +128,14 @@ int main(int argc, char *argv[])
 
         label nBadFaces = 0;
 
+        faceSet indirectPatchFaces
+        (
+            mesh,
+            "indirectPatchFaces",
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        );
+
         {
             meshFilterPtr.set(new polyMeshFilter(mesh, pointPriority));
             polyMeshFilter& meshFilter = meshFilterPtr();
@@ -139,19 +147,26 @@ int main(int argc, char *argv[])
             // the face filtering is sped up.
             nBadFaces = meshFilter.filterEdges(0);
             {
-                polyTopoChange meshMod(newMesh);
+                polyTopoChange meshMod(newMesh());
 
                 meshMod.changeMesh(mesh, false);
+
+                polyMeshFilter::copySets(newMesh(), mesh);
             }
 
             pointPriority = meshFilter.pointPriority();
         }
 
-        if (collapseFaceZone)
+        if (collapseFaceSet)
         {
-            const word faceZoneName = args.optionRead<word>("collapseFaceZone");
-
-            const faceZone& fZone = mesh.faceZones()[faceZoneName];
+            const word faceSetName(args.optionRead<word>("collapseFaceSet"));
+            faceSet fSet
+            (
+                mesh,
+                faceSetName,
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            );
 
             meshFilterPtr.reset(new polyMeshFilter(mesh, pointPriority));
             polyMeshFilter& meshFilter = meshFilterPtr();
@@ -160,11 +175,13 @@ int main(int argc, char *argv[])
 
             // Filter faces. Pass in the number of bad faces that are present
             // from the previous edge filtering to use as a stopping criterion.
-            meshFilter.filter(fZone);
+            meshFilter.filter(fSet);
             {
                 polyTopoChange meshMod(newMesh);
 
                 meshMod.changeMesh(mesh, false);
+
+                polyMeshFilter::copySets(newMesh(), mesh);
             }
 
             pointPriority = meshFilter.pointPriority();
@@ -184,6 +201,8 @@ int main(int argc, char *argv[])
                 polyTopoChange meshMod(newMesh);
 
                 meshMod.changeMesh(mesh, false);
+
+                polyMeshFilter::copySets(newMesh(), mesh);
             }
 
             pointPriority = meshFilter.pointPriority();
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/Make/options b/applications/utilities/mesh/generation/foamyHexMesh/Make/options
index f03dbf9097..f0483b99be 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/Make/options
+++ b/applications/utilities/mesh/generation/foamyHexMesh/Make/options
@@ -10,6 +10,7 @@ include $(GENERAL_RULES)/CGAL
 EXE_INC = \
     ${EXE_FROUNDING_MATH} \
     ${EXE_NDEBUG} \
+    ${CGAL_EXACT} \
     ${CGAL_INEXACT} \
     ${CGAL_INC} \
     -IconformalVoronoiMesh/lnInclude \
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
index 3afdf9cf13..8f90bb2118 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
@@ -190,20 +190,25 @@ void Foam::DelaunayMesh<Triangulation>::reset()
     resetVertexCount();
     resetCellCount();
 
-    insertPoints(vertices);
+    insertPoints(vertices, false);
 
     Info<< "Inserted " << vertexCount() << " fixed points" << endl;
 }
 
 
 template<class Triangulation>
-void Foam::DelaunayMesh<Triangulation>::insertPoints(const List<Vb>& vertices)
+Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::insertPoints
+(
+    const List<Vb>& vertices,
+    const bool reIndex
+)
 {
-    rangeInsertWithInfo
+    return rangeInsertWithInfo
     (
         vertices.begin(),
         vertices.end(),
-        false
+        false,
+        reIndex
     );
 }
 
@@ -268,7 +273,7 @@ const
 
 template<class Triangulation>
 template<class PointIterator>
-void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
+Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
 (
     PointIterator begin,
     PointIterator end,
@@ -307,6 +312,10 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
 
     Vertex_handle hint;
 
+    Map<label> oldToNewIndex(points.size());
+
+    label maxIndex = -1;
+
     for
     (
         typename vectorPairPointIndex::const_iterator p = points.begin();
@@ -335,11 +344,14 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
         {
             if (reIndex)
             {
+                const label oldIndex = vert.index();
                 hint->index() = getNewVertexIndex();
+                oldToNewIndex.insert(oldIndex, hint->index());
             }
             else
             {
                 hint->index() = vert.index();
+                maxIndex = max(maxIndex, vert.index());
             }
             hint->type() = vert.type();
             hint->procIndex() = vert.procIndex();
@@ -347,6 +359,13 @@ void Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
             hint->alignment() = vert.alignment();
         }
     }
+
+    if (!reIndex)
+    {
+        vertexCount_ = maxIndex + 1;
+    }
+
+    return oldToNewIndex;
 }
 
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H
index a7cbf26766..0bf843f1c9 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.H
@@ -228,7 +228,11 @@ public:
             void reset();
 
             //- Insert the list of vertices (calls rangeInsertWithInfo)
-            void insertPoints(const List<Vb>& vertices);
+            Map<label> insertPoints
+            (
+                const List<Vb>& vertices,
+                const bool reIndex
+            );
 
             //- Function inserting points into a triangulation and setting the
             //  index and type data of the point in the correct order. This is
@@ -237,7 +241,7 @@ public:
             //  Adapted from a post on the CGAL lists: 2010-01/msg00004.html by
             //  Sebastien Loriot (Geometry Factory).
             template<class PointIterator>
-            void rangeInsertWithInfo
+            Map<label> rangeInsertWithInfo
             (
                 PointIterator begin,
                 PointIterator end,
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
index 40c76aec9d..b598a7f5de 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DistributedDelaunayMesh.C
@@ -536,7 +536,8 @@ Foam::label Foam::DistributedDelaunayMesh<Triangulation>::referVertices
     labelPairHashSet pointsNotInserted = rangeInsertReferredWithInfo
     (
         referredVertices.begin(),
-        referredVertices.end()
+        referredVertices.end(),
+        true
     );
 
     if (!pointsNotInserted.empty())
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
index efeef35866..190b77775b 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/files
@@ -76,6 +76,7 @@ faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
 
 searchableSurfaceFeatures/searchableSurfaceFeatures.C
 searchableSurfaceFeatures/searchableBoxFeatures.C
+searchableSurfaceFeatures/searchablePlateFeatures.C
 searchableSurfaceFeatures/triSurfaceMeshFeatures.C
 
 LIB = $(FOAM_LIBBIN)/libconformalVoronoiMesh
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options
index 2be5855323..8cfdbd04fd 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/Make/options
@@ -11,6 +11,7 @@ FFLAGS = -DCGAL_FILES='"${CGAL_ARCH_PATH}/share/files"'
 EXE_INC = \
     ${EXE_FROUNDING_MATH} \
     ${EXE_NDEBUG} \
+    ${CGAL_EXACT} \
     ${CGAL_INEXACT} \
     ${CGAL_INC} \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
index 71ba6933fb..2bfe84404a 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C
@@ -628,49 +628,6 @@ Foam::tensorField Foam::cellShapeControlMesh::dumpAlignments() const
 }
 
 
-void Foam::cellShapeControlMesh::insertBoundingPoints
-(
-    const boundBox& bb,
-    const cellSizeAndAlignmentControls& sizeControls
-)
-{
-    // Loop over bound box points and get cell size and alignment
-    const pointField bbPoints(bb.points());
-
-    forAll(bbPoints, pI)
-    {
-        const Foam::point& pt = bbPoints[pI];
-
-        // Cell size here will return default cell size
-        const scalar cellSize = sizeControls.cellSize(pt);
-
-        if (debug)
-        {
-            Info<< "Insert Bounding Point: " << pt << " " << cellSize << endl;
-        }
-
-        // Get the cell size of the nearest surface.
-//        geometryToConformTo_.findSurfaceNearest
-//        (
-//            pt,
-//            GREAT,
-//            surfHit,
-//            hitSurface
-//        );
-
-        const tensor alignment = tensor::I;
-
-        insert
-        (
-            pt,
-            cellSize,
-            alignment,
-            Vb::vtInternalNearBoundary
-        );
-    }
-}
-
-
 void Foam::cellShapeControlMesh::write() const
 {
     Info<< "Writing " << meshSubDir << endl;
@@ -781,4 +738,78 @@ void Foam::cellShapeControlMesh::write() const
 }
 
 
+Foam::label Foam::cellShapeControlMesh::estimateCellCount
+(
+    const autoPtr<backgroundMeshDecomposition>& decomposition
+) const
+{
+    // Loop over all the tets and estimate the cell count in each one
+
+    scalar cellCount = 0;
+
+    for
+    (
+        Finite_cells_iterator cit = finite_cells_begin();
+        cit != finite_cells_end();
+        ++cit
+    )
+    {
+        if (!cit->hasFarPoint() && !is_infinite(cit))
+        {
+            // @todo Check if tet centre is on the processor..
+            CGAL::Tetrahedron_3<baseK> tet
+            (
+                cit->vertex(0)->point(),
+                cit->vertex(1)->point(),
+                cit->vertex(2)->point(),
+                cit->vertex(3)->point()
+            );
+
+            pointFromPoint centre = topoint(CGAL::centroid(tet));
+
+            if
+            (
+                Pstream::parRun()
+             && !decomposition().positionOnThisProcessor(centre)
+            )
+            {
+                continue;
+            }
+
+            scalar volume = CGAL::to_double(tet.volume());
+
+            scalar averagedPointCellSize = 0;
+            //scalar averagedPointCellSize = 1;
+
+            // Get an average volume by averaging the cell size of the vertices
+            for (label vI = 0; vI < 4; ++vI)
+            {
+                averagedPointCellSize += cit->vertex(vI)->targetCellSize();
+                //averagedPointCellSize *= cit->vertex(vI)->targetCellSize();
+            }
+
+            averagedPointCellSize /= 4;
+            //averagedPointCellSize = ::sqrt(averagedPointCellSize);
+
+//            if (averagedPointCellSize < SMALL)
+//            {
+//                Pout<< "Volume = " << volume << endl;
+//
+//                for (label vI = 0; vI < 4; ++vI)
+//                {
+//                    Pout<< "Point " << vI
+//                        << ", point = " << topoint(cit->vertex(vI)->point())
+//                        << ", size = " << cit->vertex(vI)->targetCellSize()
+//                        << endl;
+//                }
+//            }
+
+            cellCount += volume/pow(averagedPointCellSize, 3);
+        }
+    }
+
+    return cellCount;
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.H
index e6af5bde0f..79a397c49a 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.H
@@ -153,15 +153,14 @@ public:
 
             tensorField dumpAlignments() const;
 
-            void insertBoundingPoints
-            (
-                const boundBox& bb,
-                const cellSizeAndAlignmentControls& sizeControls
-            );
-
             void writeTriangulation();
 
             void write() const;
+
+            label estimateCellCount
+            (
+                const autoPtr<backgroundMeshDecomposition>& decomposition
+            ) const;
 };
 
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
index bff35cff69..49fc665e93 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C
@@ -429,6 +429,11 @@ void Foam::searchableSurfaceControl::initialVertices
                 vectorField normals(1);
                 searchableSurface_.getNormal(infoList, normals);
 
+                if (mag(normals[0]) < SMALL)
+                {
+                    normals[0] = vector(1, 1, 1);
+                }
+
                 pointAlignment.set(new triad(normals[0]));
 
                 if (infoList[0].hit())
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
index 8d954799a1..2be7d916d3 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/cellShapeControl/controlMeshRefinement/controlMeshRefinement.C
@@ -768,7 +768,7 @@ Foam::label Foam::controlMeshRefinement::refineMesh
         }
     }
 
-    mesh_.insertPoints(verts);
+    mesh_.insertPoints(verts, false);
 
     return verts.size();
 }
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index 10ea65f772..547a65eb12 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -179,15 +179,23 @@ void Foam::conformalVoronoiMesh::insertInternalPoints
 }
 
 
-void Foam::conformalVoronoiMesh::insertPoints
+Foam::Map<Foam::label> Foam::conformalVoronoiMesh::insertPointPairs
 (
     List<Vb>& vertices,
-    bool distribute
+    bool distribute,
+    bool reIndex
 )
 {
     if (Pstream::parRun() && distribute)
     {
-        decomposition_().distributePoints(vertices);
+        autoPtr<mapDistribute> mapDist =
+            decomposition_().distributePoints(vertices);
+
+        // Re-index the point pairs if one or both have been distributed.
+        // If both, remove
+
+        // If added a point, then need to know its point pair
+        // If one moved, then update procIndex locally
 
         forAll(vertices, vI)
         {
@@ -197,7 +205,8 @@ void Foam::conformalVoronoiMesh::insertPoints
 
     label preReinsertionSize(number_of_vertices());
 
-    this->DelaunayMesh<Delaunay>::insertPoints(vertices);
+    Map<label> oldToNewIndices =
+        this->DelaunayMesh<Delaunay>::insertPoints(vertices, reIndex);
 
     const label nReinserted = returnReduce
     (
@@ -208,17 +217,18 @@ void Foam::conformalVoronoiMesh::insertPoints
     Info<< "    Reinserted " << nReinserted << " vertices out of "
         << returnReduce(vertices.size(), sumOp<label>())
         << endl;
+
+    return oldToNewIndices;
 }
 
 
 void Foam::conformalVoronoiMesh::insertSurfacePointPairs
 (
     const pointIndexHitAndFeatureList& surfaceHits,
-    const fileName fName
+    const fileName fName,
+    DynamicList<Vb>& pts
 )
 {
-    DynamicList<Vb> pts(2.0*surfaceHits.size());
-
     forAll(surfaceHits, i)
     {
         vectorField norm(1);
@@ -246,6 +256,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
                 pointPairDistance(surfacePt),
                 surfacePt,
                 normal,
+                true,
                 pts
             );
         }
@@ -256,6 +267,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
                 pointPairDistance(surfacePt),
                 surfacePt,
                 normal,
+                true,
                 pts
             );
         }
@@ -266,6 +278,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
                 pointPairDistance(surfacePt),
                 surfacePt,
                 -normal,
+                true,
                 pts
             );
         }
@@ -280,8 +293,6 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
         }
     }
 
-    insertPoints(pts, true);
-
     if (foamyHexMeshControls().objOutput() && fName != fileName::null)
     {
         DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
@@ -292,11 +303,10 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
 void Foam::conformalVoronoiMesh::insertEdgePointGroups
 (
     const pointIndexHitAndFeatureList& edgeHits,
-    const fileName fName
+    const fileName fName,
+    DynamicList<Vb>& pts
 )
 {
-    DynamicList<Vb> pts(3.0*edgeHits.size());
-
     forAll(edgeHits, i)
     {
         if (edgeHits[i].first().hit())
@@ -318,10 +328,6 @@ void Foam::conformalVoronoiMesh::insertEdgePointGroups
         }
     }
 
-    pts.shrink();
-
-    insertPoints(pts, true);
-
     if (foamyHexMeshControls().objOutput() && fName != fileName::null)
     {
         DelaunayMeshTools::writeOBJ(time().path()/fName, pts);
@@ -859,6 +865,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
         geometryToConformTo_
     ),
     limitBounds_(),
+    ptPairs_(*this),
     ftPtConformer_(*this),
     edgeLocationTreePtr_(),
     surfacePtLocationTreePtr_(),
@@ -1108,7 +1115,7 @@ void Foam::conformalVoronoiMesh::move()
                 vB->alignment().T() & cartesianDirections
             );
 
-            Field<vector> alignmentDirs(3);
+            Field<vector> alignmentDirs(alignmentDirsA);
 
             forAll(alignmentDirsA, aA)
             {
@@ -1333,12 +1340,17 @@ void Foam::conformalVoronoiMesh::move()
                         // Point removal
                         if
                         (
-                            vA->internalPoint()
-                         && !vA->referred()
-                         && !vA->fixed()
-                         && vB->internalPoint()
-                         && !vB->referred()
-                         && !vB->fixed()
+                            (
+                                vA->internalPoint()
+                             && !vA->referred()
+                             && !vA->fixed()
+                            )
+                         &&
+                            (
+                                vB->internalPoint()
+                             && !vB->referred()
+                             && !vB->fixed()
+                            )
                         )
                         {
                             // Only insert a point at the midpoint of
@@ -1507,9 +1519,7 @@ void Foam::conformalVoronoiMesh::move()
         Info<< "Writing point displacement vectors to file." << endl;
         OFstream str
         (
-            time().path()
-          + "displacements_" + runTime_.timeName()
-          + ".obj"
+            time().path()/"displacements_" + runTime_.timeName() + ".obj"
         );
 
         for
@@ -1607,6 +1617,74 @@ void Foam::conformalVoronoiMesh::move()
                 time().path()/"boundaryPoints_" + time().timeName() + ".obj",
                 *this
             );
+
+            DelaunayMeshTools::writeOBJ
+            (
+                time().path()/"internalBoundaryPoints_" + time().timeName()
+              + ".obj",
+                *this,
+                Foam::indexedVertexEnum::vtInternalSurface,
+                Foam::indexedVertexEnum::vtInternalFeaturePoint
+            );
+
+            DelaunayMeshTools::writeOBJ
+            (
+                time().path()/"externalBoundaryPoints_" + time().timeName()
+              + ".obj",
+                *this,
+                Foam::indexedVertexEnum::vtExternalSurface,
+                Foam::indexedVertexEnum::vtExternalFeaturePoint
+            );
+
+            OBJstream multipleIntersections
+            (
+                "multipleIntersections_"
+              + time().timeName()
+              + ".obj"
+            );
+
+            for
+            (
+                Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+                eit != finite_edges_end();
+                ++eit
+            )
+            {
+                Cell_handle c = eit->first;
+                Vertex_handle vA = c->vertex(eit->second);
+                Vertex_handle vB = c->vertex(eit->third);
+
+                Foam::point ptA(topoint(vA->point()));
+                Foam::point ptB(topoint(vB->point()));
+
+                List<pointIndexHit> surfHits;
+                labelList hitSurfaces;
+
+                geometryToConformTo().findSurfaceAllIntersections
+                (
+                    ptA,
+                    ptB,
+                    surfHits,
+                    hitSurfaces
+                );
+
+                if
+                (
+                    surfHits.size() >= 2
+                 || (
+                     surfHits.size() == 0
+                  && (
+                          (vA->externalBoundaryPoint()
+                       && vB->internalBoundaryPoint())
+                       || (vB->externalBoundaryPoint()
+                       && vA->internalBoundaryPoint())
+                     )
+                    )
+                )
+                {
+                    multipleIntersections.write(linePointRef(ptA, ptB));
+                }
+            }
         }
     }
 
@@ -1629,49 +1707,6 @@ void Foam::conformalVoronoiMesh::move()
 }
 
 
-//Foam::labelListList Foam::conformalVoronoiMesh::overlapsProc
-//(
-//    const List<Foam::point>& centres,
-//    const List<scalar>& radiusSqrs
-//) const
-//{
-//    if (!Pstream::parRun())
-//    {
-//        return labelListList(centres.size(), labelList(0));
-//    }
-//
-////    DynamicList<Foam::point> pts(number_of_vertices());
-//
-////    for
-////    (
-////        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-////        vit != finite_vertices_end();
-////        vit++
-////    )
-////    {
-////        pts.append(topoint(vit->point()));
-////    }
-////
-////    dynamicIndexedOctree<dynamicTreeDataPoint> vertexOctree
-////    (
-////        dynamicTreeDataPoint(pts),
-////        treeBoundBox(min(pts), max(pts)),
-////        10, // maxLevel
-////        10, // leafSize
-////        3.0 // duplicity
-////    );
-//
-//    return decomposition_().overlapsProcessors
-//    (
-//        centres,
-//        radiusSqrs,
-//        *this,
-//        false//,
-////        vertexOctree
-//    );
-//}
-
-
 void Foam::conformalVoronoiMesh::checkCoPlanarCells() const
 {
     typedef CGAL::Exact_predicates_exact_constructions_kernel   Kexact;
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 334330f970..bc4d724ec0 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -77,7 +77,7 @@ SourceFiles
 #include "Pair.H"
 #include "DistributedDelaunayMesh.H"
 #include "featurePointConformer.H"
-
+#include "pointPairs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -86,12 +86,10 @@ namespace Foam
 
 // Forward declaration of classes
 class initialPointsMethod;
-
 class relaxationModel;
-
 class faceAreaWeightModel;
-
 class backgroundMeshDecomposition;
+class OBJstream;
 
 /*---------------------------------------------------------------------------*\
                     Class conformalVoronoiMesh Declaration
@@ -165,6 +163,9 @@ private:
         //- Limiting bound box before infinity begins
         treeBoundBox limitBounds_;
 
+        //-
+        mutable pointPairs<Delaunay> ptPairs_;
+
         //-
         featurePointConformer ftPtConformer_;
 
@@ -229,10 +230,11 @@ private:
             const bool distribute = false
         );
 
-        void insertPoints
+        Map<label> insertPointPairs
         (
             List<Vb>& vertices,
-            bool distribute
+            bool distribute,
+            bool reIndex
         );
 
         //- Create a point-pair at a ppDist distance either side of
@@ -242,6 +244,7 @@ private:
             const scalar ppDist,
             const Foam::point& surfPt,
             const vector& n,
+            const bool ptPair,
             DynamicList<Vb>& pts
         ) const;
 
@@ -254,24 +257,20 @@ private:
             const scalar ppDist,
             const Foam::point& surfPt,
             const vector& n,
+            const bool ptPair,
             DynamicList<Vb>& pts
         ) const;
 
         //- Check internal point is completely inside the meshable region
         inline bool internalPointIsInside(const Foam::point& pt) const;
 
-        inline bool isPointPair
-        (
-            const Vertex_handle& vA,
-            const Vertex_handle& vB
-        ) const;
-
         //- Insert pairs of points on the surface with the given normals, at the
         //  specified spacing
         void insertSurfacePointPairs
         (
             const pointIndexHitAndFeatureList& surfaceHits,
-            const fileName fName = fileName::null
+            const fileName fName,
+            DynamicList<Vb>& pts
         );
 
         //- Insert groups of points to conform to an edge given a list of
@@ -280,7 +279,8 @@ private:
         void insertEdgePointGroups
         (
             const pointIndexHitAndFeatureList& edgeHits,
-            const fileName fName = fileName::null
+            const fileName fName,
+            DynamicList<Vb>& pts
         );
 
         void createEdgePointGroupByCirculating
@@ -474,14 +474,11 @@ private:
             label& hitSurface
         ) const;
 
-        //- Find the "worst" incursion of the dual cell of a non-internal or
-        //  boundary point through the surface, subject to the
-        //  maxSurfaceProtrusion tolerance
         void dualCellLargestSurfaceIncursion
         (
             const Delaunay::Finite_vertices_iterator& vit,
-            pointIndexHit& surfHitLargest,
-            label& hitSurfaceLargest
+            pointIndexHit& surfHit,
+            label& hitSurface
         ) const;
 
         //- Write out vertex-processor occupancy information for debugging
@@ -695,7 +692,8 @@ private:
         //- Merge vertices that are identical
         void mergeIdenticalDualVertices
         (
-            const pointField& pts
+            const pointField& pts,
+            labelList& boundaryPts
         );
 
         label mergeIdenticalDualVertices
@@ -736,7 +734,11 @@ private:
         );
 
         //- Re-index all of the Delaunay cells
-        void reindexDualVertices(const Map<label>& dualPtIndexMap);
+        void reindexDualVertices
+        (
+            const Map<label>& dualPtIndexMap,
+            labelList& boundaryPts
+        );
 
         label createPatchInfo
         (
@@ -846,6 +848,8 @@ private:
             const PtrList<dictionary>& patchDicts
         ) const;
 
+        void writePointPairs(const fileName& fName) const;
+
         //- Disallow default bitwise copy construct
         conformalVoronoiMesh(const conformalVoronoiMesh&);
 
@@ -992,7 +996,7 @@ public:
                 const wordList& patchNames,
                 const PtrList<dictionary>& patchDicts,
                 const pointField& cellCentres,
-                const PackedBoolList& boundaryFacesToRemove
+                PackedBoolList& boundaryFacesToRemove
             ) const;
 
             //- Calculate and write a field of the target cell size,
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index cb17a1367d..e220556fa5 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -30,487 +30,11 @@ License
 #include "indexedCellChecks.H"
 #include "OBJstream.H"
 #include "indexedCellOps.H"
+#include "ListOps.H"
 #include "DelaunayMeshTools.H"
 
-#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
-#include "CGAL/Gmpq.h"
-
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-void Foam::conformalVoronoiMesh::checkCells()
-{
-    List<List<FixedList<Foam::point, 4> > > cellListList(Pstream::nProcs());
-
-    List<FixedList<Foam::point, 4> > cells(number_of_finite_cells());
-
-    globalIndex gIndex(number_of_vertices());
-
-    label count = 0;
-    for
-    (
-        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-        cit != finite_cells_end();
-        ++cit
-    )
-    {
-        if (tetrahedron(cit).volume() == 0)
-        {
-            Pout<< "ZERO VOLUME TET" << endl;
-            Pout<< cit->info();
-            Pout<< cit->dual();
-        }
-
-        if (cit->hasFarPoint())
-        {
-            continue;
-        }
-
-        List<labelPair> cellVerticesPair(4);
-        List<Foam::point> cellVertices(4);
-
-        for (label vI = 0; vI < 4; ++vI)
-        {
-            cellVerticesPair[vI] = labelPair
-            (
-                cit->vertex(vI)->procIndex(),
-                cit->vertex(vI)->index()
-            );
-            cellVertices[vI] = topoint(cit->vertex(vI)->point());
-        }
-
-        List<Foam::point> cellVerticesOld(cellVertices);
-        labelList oldToNew;
-        sortedOrder(cellVerticesPair, oldToNew);
-        oldToNew = invert(oldToNew.size(), oldToNew);
-        inplaceReorder(oldToNew, cellVerticesPair);
-        inplaceReorder(oldToNew, cellVertices);
-
-//        FixedList<label, 4> globalTetCell
-//        (
-//            cit->globallyOrderedCellVertices(gIndex)
-//        );
-//
-//        FixedList<Point, 4> cellVertices(Point(0,0,0));
-//
-//        forAll(globalTetCell, gvI)
-//        {
-//            label gI = globalTetCell[gvI];
-//
-//            cellVertices[gvI] = cit->vertex(gI)->point();
-//        }
-
-//        if (cit->hasFarPoint())
-//        {
-//            continue;
-//        }
-
-        for (label i = 0; i < 4; ++i)
-        {
-            //cells[count][i] = topoint(cit->vertex(i)->point());
-            cells[count][i] = cellVertices[i];
-        }
-
-        count++;
-    }
-
-    cells.setSize(count);
-
-    cellListList[Pstream::myProcNo()] = cells;
-
-    Pstream::gatherList(cellListList);
-
-    if (Pstream::master())
-    {
-        Info<< "Checking on master processor the cells of each " << nl
-            << "processor point list against the master cell list." << nl
-            << "There are " << cellListList.size() << " processors" << nl
-            << "The size of each processor's cell list is:" << endl;
-
-        forAll(cellListList, cfI)
-        {
-            Info<< "    Proc " << cfI << " has " << cellListList[cfI].size()
-                << " cells" << endl;
-        }
-
-        label nMatches = 0, nMatchFoundDiffOrder = 0;
-
-        forAll(cellListList[0], cmI)
-        {
-            const FixedList<Foam::point, 4>& masterCell = cellListList[0][cmI];
-
-            bool matchFound = false;
-            bool matchFoundDiffOrder = false;
-
-            forAll(cellListList, cpI)
-            {
-                if (cpI == 0)
-                {
-                    continue;
-                }
-
-                forAll(cellListList[cpI], csI)
-                {
-                    const FixedList<Foam::point, 4>& slaveCell
-                        = cellListList[cpI][csI];
-
-                    if (masterCell == slaveCell)
-                    {
-                        matchFound = true;
-                        break;
-                    }
-                    else
-                    {
-                        label samePt = 0;
-
-                        forAll(masterCell, mI)
-                        {
-                            const Foam::point& mPt = masterCell[mI];
-
-                            forAll(slaveCell, sI)
-                            {
-                                const Foam::point& sPt = slaveCell[sI];
-
-                                if (mPt == sPt)
-                                {
-                                    samePt++;
-                                }
-                            }
-                        }
-
-                        if (samePt == 4)
-                        {
-                            matchFoundDiffOrder = true;
-
-                            Pout<< masterCell << nl << slaveCell << endl;
-
-                            break;
-                        }
-                    }
-                }
-            }
-
-            if (matchFound)
-            {
-                nMatches++;
-            }
-
-            if (matchFoundDiffOrder)
-            {
-                nMatchFoundDiffOrder++;
-            }
-        }
-
-        Info<< "Found " << nMatches << " matching cells and "
-            << nMatchFoundDiffOrder << " matching cells with different "
-            << "vertex ordering"<< endl;
-    }
-}
-
-
-//void Foam::conformalVoronoiMesh::checkDuals()
-//{
-//    List<List<Point> > pointFieldList(Pstream::nProcs());
-//
-//    List<Point> duals(number_of_finite_cells());
-//
-//    typedef CGAL::Exact_predicates_exact_constructions_kernel       EK2;
-//    typedef CGAL::Regular_triangulation_euclidean_traits_3<EK2>     EK;
-//    typedef CGAL::Cartesian_converter<baseK::Kernel, EK2>  To_exact;
-//    typedef CGAL::Cartesian_converter<EK2, baseK::Kernel>  Back_from_exact;
-//
-////    PackedBoolList bPoints(number_of_finite_cells());
-//
-////    indexDualVertices(duals, bPoints);
-//
-//    label count = 0;//duals.size();
-//
-//    duals.setSize(number_of_finite_cells());
-//
-//    globalIndex gIndex(number_of_vertices());
-//
-//    for
-//    (
-//        Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-//        cit != finite_cells_end();
-//        ++cit
-//    )
-//    {
-//        if (cit->hasFarPoint())
-//        {
-//            continue;
-//        }
-//
-//        duals[count++] = cit->circumcenter();
-//
-////        List<labelPair> cellVerticesPair(4);
-////        List<Point> cellVertices(4);
-////
-////        for (label vI = 0; vI < 4; ++vI)
-////        {
-////            cellVerticesPair[vI] = labelPair
-////            (
-////                cit->vertex(vI)->procIndex(),
-////                cit->vertex(vI)->index()
-////            );
-////            cellVertices[vI] = cit->vertex(vI)->point();
-////        }
-////
-////        labelList oldToNew;
-////        sortedOrder(cellVerticesPair, oldToNew);
-////        oldToNew = invert(oldToNew.size(), oldToNew);
-////        inplaceReorder(oldToNew, cellVerticesPair);
-////        inplaceReorder(oldToNew, cellVertices);
-////
-////        duals[count++] = CGAL::circumcenter
-////        (
-////            cellVertices[0],
-////            cellVertices[1],
-////            cellVertices[2],
-////            cellVertices[3]
-////        );
-//
-////        To_exact to_exact;
-////        Back_from_exact back_from_exact;
-////        EK::Construct_circumcenter_3 exact_circumcenter =
-////            EK().construct_circumcenter_3_object();
-////
-////        duals[count++] = topoint
-////        (
-////            back_from_exact
-////            (
-////                exact_circumcenter
-////                (
-////                    to_exact(cit->vertex(0)->point()),
-////                    to_exact(cit->vertex(1)->point()),
-////                    to_exact(cit->vertex(2)->point()),
-////                    to_exact(cit->vertex(3)->point())
-////                )
-////            )
-////        );
-//    }
-//
-//    Pout<< "Duals Calculated " << count << endl;
-//
-//    duals.setSize(count);
-//
-//    pointFieldList[Pstream::myProcNo()] = duals;
-//
-//    Pstream::gatherList(pointFieldList);
-//
-//    if (Pstream::master())
-//    {
-//        Info<< "Checking on master processor the dual locations of each" << nl
-//            << " processor point list against the master dual list." << nl
-//            << "There are " << pointFieldList.size() << " processors" << nl
-//            << "The size of each processor's dual list is:" << endl;
-//
-//        forAll(pointFieldList, pfI)
-//        {
-//            Info<< "    Proc " << pfI << " has " << pointFieldList[pfI].size()
-//                << " duals" << endl;
-//        }
-//
-//        label nNonMatches = 0;
-//        label nNearMatches = 0;
-//        label nExactMatches = 0;
-//
-//        forAll(pointFieldList[0], pI)
-//        {
-//            const Point& masterPoint = pointFieldList[0][pI];
-//
-//            bool foundMatch = false;
-//            bool foundNearMatch = false;
-//
-//            scalar minCloseness = GREAT;
-//            Point closestPoint(0, 0, 0);
-//
-//            forAll(pointFieldList, pfI)
-//            {
-//                if (pfI == 0)
-//                {
-//                    continue;
-//                }
-//
-////                label pfI = 1;
-//
-//                forAll(pointFieldList[pfI], pISlave)
-//                {
-//                    const Point& slavePoint
-//                        = pointFieldList[pfI][pISlave];
-//
-//                    if (masterPoint == slavePoint)
-//                    {
-//                        foundMatch = true;
-//                        break;
-//                    }
-//
-//                    const scalar closeness = mag
-//                    (
-//                        topoint(masterPoint) - topoint(slavePoint)
-//                    );
-//
-//                    if (closeness < 1e-12)
-//                    {
-//                        foundNearMatch = true;
-//                    }
-//                    else
-//                    {
-//                        if (closeness < minCloseness)
-//                        {
-//                            minCloseness = closeness;
-//                            closestPoint = slavePoint;
-//                        }
-//                    }
-//                }
-//
-//                if (!foundMatch)
-//                {
-//                    if (foundNearMatch)
-//                    {
-//                        CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
-//                        CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
-//                        CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
-//
-//                        std::cout<< "master = " << x << " " << y << " " << z
-//                            << std::endl;
-//
-//                        CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
-//                        CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
-//                        CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
-//                        std::cout<< "slave  = " << xs << " " << ys << " "
-//                                 << zs
-//                            << std::endl;
-//
-//                        nNearMatches++;
-//                    }
-//                    else
-//                    {
-//                        nNonMatches++;
-//                        Info<< "Closest point to " << masterPoint << " is "
-//                            << closestPoint << nl
-//                            << "    Separation is " << minCloseness << endl;
-//
-//                        CGAL::Gmpq x(CGAL::to_double(masterPoint.x()));
-//                        CGAL::Gmpq y(CGAL::to_double(masterPoint.y()));
-//                        CGAL::Gmpq z(CGAL::to_double(masterPoint.z()));
-//
-//                        std::cout<< "master = " << x << " " << y << " " << z
-//                                 << std::endl;
-//
-//                        CGAL::Gmpq xs(CGAL::to_double(closestPoint.x()));
-//                        CGAL::Gmpq ys(CGAL::to_double(closestPoint.y()));
-//                        CGAL::Gmpq zs(CGAL::to_double(closestPoint.z()));
-//                        std::cout<< "slave  = " << xs << " " << ys << " "
-//                                 << zs
-//                                 << std::endl;
-//                    }
-//                }
-//                else
-//                {
-//                    nExactMatches++;
-//                }
-//            }
-//        }
-//
-//        Info<< "Found " << nNonMatches << " non-matching duals" << nl
-//            << " and " << nNearMatches << " near matches"
-//            << " and " << nExactMatches << " exact matches" << endl;
-//    }
-//}
-
-
-void Foam::conformalVoronoiMesh::checkVertices()
-{
-    List<pointField> pointFieldList(Pstream::nProcs());
-
-    pointField points(number_of_vertices());
-
-    labelPairHashSet duplicateVertices;
-
-    label count = 0;
-    for
-    (
-        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-        vit != finite_vertices_end();
-        ++vit
-    )
-    {
-        if (duplicateVertices.found(labelPair(vit->procIndex(), vit->index())))
-        {
-            Pout<< "DUPLICATE " << vit->procIndex() << vit->index() << endl;
-        }
-        else
-        {
-            duplicateVertices.insert(labelPair(vit->procIndex(), vit->index()));
-        }
-
-        points[count++] = topoint(vit->point());
-    }
-
-    pointFieldList[Pstream::myProcNo()] = points;
-
-    Pstream::gatherList(pointFieldList);
-
-    OFstream str("missingPoints.obj");
-
-    if (Pstream::master())
-    {
-        Info<< "Checking on master processor the point locations of each " << nl
-            << "processor point list against the master point list." << nl
-            << "There are " << pointFieldList.size() << " processors" << nl
-            << "The size of each processor's point list is:" << endl;
-
-        forAll(pointFieldList, pfI)
-        {
-            Info<< "    Proc " << pfI << " has " << pointFieldList[pfI].size()
-                << " points" << endl;
-        }
-
-        label nNonMatches = 0;
-
-        forAll(pointFieldList[0], pI)
-        {
-            const Foam::point& masterPoint = pointFieldList[0][pI];
-
-            forAll(pointFieldList, pfI)
-            {
-                if (pI == 0)
-                {
-                    continue;
-                }
-
-                bool foundMatch = false;
-
-                forAll(pointFieldList[pfI], pISlave)
-                {
-                    const Foam::point& slavePoint
-                        = pointFieldList[pfI][pISlave];
-
-                    if (masterPoint == slavePoint)
-                    {
-                        foundMatch = true;
-                        break;
-                    }
-                }
-
-                if (!foundMatch)
-                {
-                    Info<< "    Proc " << pfI << " Master != Slave -> "
-                        << masterPoint << endl;
-
-                    meshTools::writeOBJ(str, masterPoint);
-
-                    nNonMatches++;
-                }
-            }
-        }
-
-        Info<< "Found a total of " << nNonMatches << " non-matching points"
-            << endl;
-    }
-}
-
-
 void Foam::conformalVoronoiMesh::calcDualMesh
 (
     pointField& points,
@@ -528,53 +52,6 @@ void Foam::conformalVoronoiMesh::calcDualMesh
 {
     timeCheck("Start calcDualMesh");
 
-//    if (debug)
-//    {
-//        Pout<< nl << "Perfoming some checks . . ." << nl << nl
-//            << "Total number of vertices = " << number_of_vertices() << nl
-//            << "Total number of cells    = " << number_of_finite_cells()
-//            << endl;
-//
-//        checkVertices();
-//        checkCells();
-//        checkDuals();
-//
-//        Info<< nl << "Finished checks" << nl << endl;
-//    }
-
-//    OFstream str("attachedToFeature.obj");
-//    label offset = 0;
-//
-//    for
-//    (
-//        Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
-//        vit != finite_vertices_end();
-//        ++vit
-//    )
-//    {
-//        if (vit->featurePoint())
-//        {
-//            std::list<Cell_handle> adjacentCells;
-//
-//            finite_incident_cells(vit, std::back_inserter(adjacentCells));
-//
-//            for
-//            (
-//                std::list<Cell_handle>::iterator acit = adjacentCells.begin();
-//                acit != adjacentCells.end();
-//                ++acit
-//            )
-//            {
-//                if ((*acit)->real())
-//                {
-//                    drawDelaunayCell(str, (*acit), offset);
-//                    offset++;
-////                    meshTools::writeOBJ(str, topoint((*acit)->dual()));
-//                }
-//            }
-//        }
-//    }
-
     setVertexSizeAndAlignment();
 
     timeCheck("After setVertexSizeAndAlignment");
@@ -585,7 +62,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
         Info<< nl << "Merging identical points" << endl;
 
         // There is no guarantee that a merge of close points is no-risk
-        mergeIdenticalDualVertices(points);
+        mergeIdenticalDualVertices(points, boundaryPts);
     }
 
     // Final dual face and owner neighbour construction
@@ -820,7 +297,8 @@ void Foam::conformalVoronoiMesh::calcTetMesh
 
 void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
 (
-    const pointField& pts
+    const pointField& pts,
+    labelList& boundaryPts
 )
 {
     // Assess close points to be merged
@@ -838,7 +316,7 @@ void Foam::conformalVoronoiMesh::mergeIdenticalDualVertices
             dualPtIndexMap
         );
 
-        reindexDualVertices(dualPtIndexMap);
+        reindexDualVertices(dualPtIndexMap, boundaryPts);
 
         reduce(nPtsMerged, sumOp<label>());
 
@@ -1245,7 +723,6 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
         false
     );
 
-    //createCellCentres(cellCentres);
     cellCentres = DelaunayMeshTools::allPoints(*this);
 
     labelList cellToDelaunayVertex(removeUnusedCells(owner, neighbour));
@@ -1327,24 +804,6 @@ Foam::conformalVoronoiMesh::createPolyMeshFromPoints
 
     pMesh.addPatches(patches);
 
-    // Info<< "ADDPATCHES NOT IN PARALLEL" << endl;
-
-    // forAll(patches, p)
-    // {
-    //     patches[p] = new polyPatch
-    //     (
-    //         patchNames[p],
-    //         patchSizes[p],
-    //         patchStarts[p],
-    //         p,
-    //         pMesh.boundaryMesh()
-    //     );
-    // }
-
-    // pMesh.addPatches(patches, false);
-
-    // pMesh.overrideCellCentres(cellCentres);
-
     return meshPtr;
 }
 
@@ -1361,7 +820,7 @@ void Foam::conformalVoronoiMesh::checkCellSizing()
     indexDualVertices(ptsField, boundaryPts);
 
     // Merge close dual vertices.
-    mergeIdenticalDualVertices(ptsField);
+    mergeIdenticalDualVertices(ptsField, boundaryPts);
 
     autoPtr<polyMesh> meshPtr = createPolyMeshFromPoints(ptsField);
     const polyMesh& pMesh = meshPtr();
@@ -2004,23 +1463,10 @@ void Foam::conformalVoronoiMesh::indexDualVertices
                     boundaryPts[cit->cellIndex()] = surface;
                 }
             }
-            else if
-            (
-                cit->baffleBoundaryDualVertex()
-            )
+            else if (cit->baffleBoundaryDualVertex())
             {
                 boundaryPts[cit->cellIndex()] = surface;
             }
-            else if
-            (
-                cit->vertex(0)->featureEdgePoint()
-             && cit->vertex(1)->featureEdgePoint()
-             && cit->vertex(2)->featureEdgePoint()
-             && cit->vertex(3)->featureEdgePoint()
-            )
-            {
-                boundaryPts[cit->cellIndex()] = featureEdge;
-            }
             else
             {
                 boundaryPts[cit->cellIndex()] = internal;
@@ -2040,7 +1486,8 @@ void Foam::conformalVoronoiMesh::indexDualVertices
 
 void Foam::conformalVoronoiMesh::reindexDualVertices
 (
-    const Map<label>& dualPtIndexMap
+    const Map<label>& dualPtIndexMap,
+    labelList& boundaryPts
 )
 {
     for
@@ -2053,6 +1500,12 @@ void Foam::conformalVoronoiMesh::reindexDualVertices
         if (dualPtIndexMap.found(cit->cellIndex()))
         {
             cit->cellIndex() = dualPtIndexMap[cit->cellIndex()];
+            boundaryPts[cit->cellIndex()] =
+                max
+                (
+                    boundaryPts[cit->cellIndex()],
+                    boundaryPts[dualPtIndexMap[cit->cellIndex()]]
+                );
         }
     }
 }
@@ -2267,11 +1720,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
     bool includeEmptyPatches
 ) const
 {
-    const label defaultPatchIndex = createPatchInfo
-    (
-        patchNames,
-        patchDicts
-    );
+    const label defaultPatchIndex = createPatchInfo(patchNames, patchDicts);
 
     const label nPatches = patchNames.size();
 
@@ -2296,6 +1745,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 
     List<DynamicList<bool> > indirectPatchFace(nPatches, DynamicList<bool>(0));
 
+
     faces.setSize(number_of_finite_edges());
     owner.setSize(number_of_finite_edges());
     neighbour.setSize(number_of_finite_edges());
@@ -2765,7 +2215,13 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 
                     // If the two vertices are a pair, then the patch face is
                     // a desired one.
-                    if (!isPointPair(vA, vB))
+                    if
+                    (
+                        vA->boundaryPoint() && vB->boundaryPoint()
+                     && !isProcBoundaryEdge(eit)
+                     && !ptPairs_.isPointPair(vA, vB)
+                     && !ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
+                    )
                     {
                         indirectPatchFace[patchIndex].append(true);
                     }
@@ -2797,6 +2253,18 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 
                     if (patchIndex != -1)
                     {
+//                        if
+//                        (
+//                            vA->boundaryPoint() && vB->boundaryPoint()
+//                         && !ptPairs_.isPointPair(vA, vB)
+//                        )
+//                        {
+//                            indirectPatchFace[patchIndex].append(true);
+//                        }
+//                        else
+//                        {
+//                            indirectPatchFace[patchIndex].append(false);
+//                        }
 //                        patchFaces[patchIndex].append(newDualFace);
 //                        patchOwners[patchIndex].append(own);
 //                        indirectPatchFace[patchIndex].append(false);
@@ -2813,8 +2281,6 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
 //                        {
 //                            patchPPSlaves[patchIndex].append(vA->index());
 //                        }
-
-//                        baffleFaces[dualFaceI] = patchIndex;
                     }
 //                    else
                     {
@@ -2989,7 +2455,6 @@ void Foam::conformalVoronoiMesh::sortProcPatches
         faceList& faces = patchFaces[patchI];
         labelList& owner = patchOwners[patchI];
         DynamicList<label>& slaves = patchPointPairSlaves[patchI];
-
         DynamicList<Pair<labelPair> >& sortingIndices
             = patchSortingIndices[patchI];
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index f5a023592c..b666539262 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -62,11 +62,13 @@ void Foam::conformalVoronoiMesh::conformToSurface()
 
         if (Pstream::parRun())
         {
-            sync(decomposition_().procBounds());
+            sync(decomposition().procBounds());
         }
     }
     else
     {
+        ptPairs_.clear();
+
         // Rebuild, insert and store new surface conformation
         buildSurfaceConformation();
 
@@ -74,7 +76,7 @@ void Foam::conformalVoronoiMesh::conformToSurface()
         {
             if (Pstream::parRun())
             {
-                sync(decomposition_().procBounds());
+                sync(decomposition().procBounds());
             }
         }
 
@@ -358,10 +360,13 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
         }
 
+        DynamicList<Vb> pts(2*surfaceHits.size() + 3*featureEdgeHits.size());
+
         insertSurfacePointPairs
         (
             surfaceHits,
-            "surfaceConformationLocations_initial.obj"
+            "surfaceConformationLocations_initial.obj",
+            pts
         );
 
         // In parallel, synchronise the edge trees
@@ -373,9 +378,21 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
         insertEdgePointGroups
         (
             featureEdgeHits,
-            "edgeConformationLocations_initial.obj"
+            "edgeConformationLocations_initial.obj",
+            pts
         );
 
+        pts.shrink();
+
+        Map<label> oldToNewIndices = insertPointPairs(pts, true, true);
+
+        // Re-index the point pairs
+        ptPairs_.reIndex(oldToNewIndices);
+
+        writePointPairs("pointPairs_initial.obj");
+
+        // Remove location from surface/edge tree
+
         timeCheck("After initial conformation");
 
         initialTotalHits = nSurfHits + nFeatEdHits;
@@ -550,215 +567,6 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             }
         }
 
-//        for
-//        (
-//            Delaunay::Finite_edges_iterator eit = finite_edges_begin();
-//            eit != finite_edges_end();
-//            ++eit
-//        )
-//        {
-//            Cell_handle c = eit->first;
-//            Vertex_handle vA = c->vertex(eit->second);
-//            Vertex_handle vB = c->vertex(eit->third);
-//
-//            if
-//            (
-//                vA->referred()
-//             || vB->referred()
-//            )
-//            {
-//                continue;
-//            }
-//
-//            if
-//            (
-//                (vA->internalPoint() && vA->referred())
-//             || (vB->internalPoint() && vB->referred())
-//            )
-//            {
-//                continue;
-//            }
-//
-//            if
-//            (
-//                (vA->internalPoint() && vB->externalBoundaryPoint())
-//             || (vB->internalPoint() && vA->externalBoundaryPoint())
-//             || (vA->internalBoundaryPoint() && vB->internalBoundaryPoint())
-//            )
-//            {
-//                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
-//                pointIndexHit surfHit;
-//                label hitSurface;
-//
-//                geometryToConformTo_.findSurfaceNearest
-//                (
-//                    (
-//                        vA->internalPoint()
-//                      ? topoint(vA->point())
-//                      : topoint(vB->point())
-//                    ),
-//                    magSqr(topoint(vA->point()) - topoint(vB->point())),
-//                    surfHit,
-//                    hitSurface
-//                );
-//
-//                if (surfHit.hit())
-//                {
-//                    surfaceIntersections.append
-//                    (
-//                        pointIndexHitAndFeature(surfHit, hitSurface)
-//                    );
-//
-//                    addSurfaceAndEdgeHits
-//                    (
-//                        (
-//                            vA->internalPoint()
-//                          ? topoint(vA->point())
-//                          : topoint(vB->point())
-//                        ),
-//                        surfaceIntersections,
-//                        surfacePtReplaceDistCoeffSqr,
-//                        edgeSearchDistCoeffSqr,
-//                        surfaceHits,
-//                        featureEdgeHits,
-//                        surfaceToTreeShape,
-//                        edgeToTreeShape,
-//                        surfacePtToEdgePtDist,
-//                        false
-//                    );
-//                }
-//            }
-//        }
-
-        for
-        (
-            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-            cit != finite_cells_end();
-            ++cit
-        )
-        {
-            label nInternal = 0;
-            for (label vI = 0; vI < 4; vI++)
-            {
-                if (cit->vertex(vI)->internalPoint())
-                {
-                    nInternal++;
-                }
-            }
-
-            if (nInternal == 1 && cit->hasBoundaryPoint())
-            //if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
-            {
-                const Foam::point& pt = cit->dual();
-
-                const Foam::point cellCentre =
-                    topoint
-                    (
-                        CGAL::centroid
-                        (
-                            CGAL::Tetrahedron_3<baseK>
-                            (
-                                cit->vertex(0)->point(),
-                                cit->vertex(1)->point(),
-                                cit->vertex(2)->point(),
-                                cit->vertex(3)->point()
-                            )
-                        )
-                    );
-
-                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
-                pointIndexHit surfHit;
-                label hitSurface;
-
-                geometryToConformTo_.findSurfaceNearestIntersection
-                (
-                    cellCentre,
-                    pt,
-                    surfHit,
-                    hitSurface
-                );
-
-                if (surfHit.hit())
-                {
-                    surfaceIntersections.append
-                    (
-                        pointIndexHitAndFeature(surfHit, hitSurface)
-                    );
-
-                    addSurfaceAndEdgeHits
-                    (
-                        pt,
-                        surfaceIntersections,
-                        surfacePtReplaceDistCoeffSqr,
-                        edgeSearchDistCoeffSqr,
-                        surfaceHits,
-                        featureEdgeHits,
-                        surfaceToTreeShape,
-                        edgeToTreeShape,
-                        surfacePtToEdgePtDist,
-                        false
-                    );
-                }
-            }
-        }
-
-//        for
-//        (
-//            Delaunay::Finite_cells_iterator cit = finite_cells_begin();
-//            cit != finite_cells_end();
-//            ++cit
-//        )
-//        {
-//            if (cit->boundaryDualVertex() && !cit->hasReferredPoint())
-//            {
-//                const Foam::point& pt = cit->dual();
-//
-//                pointIndexHitAndFeatureDynList surfaceIntersections(0.5*AtoV);
-//                pointIndexHit surfHit;
-//                label hitSurface;
-//
-//                geometryToConformTo_.findSurfaceNearest
-//                (
-//                    pt,
-//                    sqr(0.1*targetCellSize(pt)),
-//                    surfHit,
-//                    hitSurface
-//                );
-//
-//                if (!surfHit.hit())
-//                {
-//                    geometryToConformTo_.findSurfaceNearest
-//                    (
-//                        pt,
-//                        GREAT,
-//                        surfHit,
-//                        hitSurface
-//                    );
-//
-//                    if (surfHit.hit())
-//                    {
-//                        surfaceIntersections.append
-//                        (
-//                            pointIndexHitAndFeature(surfHit, hitSurface)
-//                        );
-//
-//                        addSurfaceAndEdgeHits
-//                        (
-//                            pt,
-//                            surfaceIntersections,
-//                            surfacePtReplaceDistCoeffSqr,
-//                            edgeSearchDistCoeffSqr,
-//                            surfaceHits,
-//                            featureEdgeHits,
-//                            surfaceToTreeShape,
-//                            edgeToTreeShape,
-//                            false
-//                        );
-//                    }
-//                }
-//            }
-//        }
-
         label nVerts = number_of_vertices();
         label nSurfHits = surfaceHits.size();
         label nFeatEdHits = featureEdgeHits.size();
@@ -789,10 +597,16 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
                     synchroniseSurfaceTrees(surfaceToTreeShape, surfaceHits);
             }
 
+            DynamicList<Vb> pts
+            (
+                2*surfaceHits.size() + 3*featureEdgeHits.size()
+            );
+
             insertSurfacePointPairs
             (
                 surfaceHits,
-                "surfaceConformationLocations_" + name(iterationNo) + ".obj"
+                "surfaceConformationLocations_" + name(iterationNo) + ".obj",
+                pts
             );
 
             // In parallel, synchronise the edge trees
@@ -805,9 +619,19 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation()
             insertEdgePointGroups
             (
                 featureEdgeHits,
-                "edgeConformationLocations_" + name(iterationNo) + ".obj"
+                "edgeConformationLocations_" + name(iterationNo) + ".obj",
+                pts
             );
 
+            pts.shrink();
+
+            Map<label> oldToNewIndices = insertPointPairs(pts, true, true);
+
+            // Reindex the point pairs
+            ptPairs_.reIndex(oldToNewIndices);
+
+            writePointPairs("pointPairs_" + name(iterationNo) + ".obj");
+
             if (Pstream::parRun())
             {
                 sync
@@ -1155,8 +979,6 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
          || is_infinite(fit->first->neighbor(fit->second))
          || !fit->first->hasInternalPoint()
          || !fit->first->neighbor(fit->second)->hasInternalPoint()
-//         || fit->first->hasFarPoint()
-//         || fit->first->neighbor(fit->second)->hasFarPoint()
         )
         {
             continue;
@@ -1209,14 +1031,13 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
 
             const scalar d = p.normalIntersect(r);
 
-            const Foam::point newPoint = vertex + d*n;
+            Foam::point newPoint = vertex + d*n;
 
             pointIndexHitAndFeature info;
-
             geometryToConformTo_.findSurfaceNearest
             (
                 newPoint,
-                surfaceSearchDistanceSqr(newPoint),
+                4.0*magSqr(newPoint - vertex),
                 info.first(),
                 info.second()
             );
@@ -1359,7 +1180,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
         if
         (
             is_infinite(c1) || is_infinite(c2)
-         || (!c1->hasInternalPoint() && !c2->hasInternalPoint())
+         || (
+                !c1->internalOrBoundaryDualVertex()
+             || !c2->internalOrBoundaryDualVertex()
+            )
          || !c1->real() || !c2->real()
         )
         {
@@ -1369,19 +1193,24 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
 //        Foam::point endPt = 0.5*(c1->dual() + c2->dual());
         Foam::point endPt = c1->dual();
 
-        if
-        (
-            magSqr(vert - c1->dual())
-          < magSqr(vert - c2->dual())
-        )
+        if (magSqr(vert - c1->dual()) < magSqr(vert - c2->dual()))
         {
             endPt = c2->dual();
         }
 
+        if
+        (
+            magSqr(vert - endPt)
+          > magSqr(geometryToConformTo().globalBounds().mag())
+        )
+        {
+            continue;
+        }
+
         pointIndexHit surfHit;
         label hitSurface;
 
-        geometryToConformTo_.findSurfaceAnyIntersection
+        geometryToConformTo_.findSurfaceNearestIntersection
         (
             vert,
             endPt,
@@ -1401,15 +1230,46 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
 
             const vector& n = norm[0];
 
-            const scalar normalProtrusionDistance =
-                (endPt - surfHit.hitPoint()) & n;
+            const scalar normalProtrusionDistance
+            (
+                (endPt - surfHit.hitPoint()) & n
+            );
 
             if (normalProtrusionDistance > maxProtrusionDistance)
             {
-                surfHitLargest = surfHit;
-                hitSurfaceLargest = hitSurface;
+                const plane p(surfHit.hitPoint(), n);
 
-                maxProtrusionDistance = normalProtrusionDistance;
+                const plane::ray r(endPt, -n);
+
+                const scalar d = p.normalIntersect(r);
+
+                Foam::point newPoint = endPt - d*n;
+
+                pointIndexHitAndFeature info;
+                geometryToConformTo_.findSurfaceNearest
+                (
+                    newPoint,
+                    4.0*magSqr(newPoint - endPt),
+                    info.first(),
+                    info.second()
+                );
+
+                if (info.first().hit())
+                {
+                    if
+                    (
+                        surfaceLocationConformsToInside
+                        (
+                            pointIndexHitAndFeature(info.first(), info.second())
+                        )
+                    )
+                    {
+                        surfHitLargest = info.first();
+                        hitSurfaceLargest = info.second();
+
+                        maxProtrusionDistance = normalProtrusionDistance;
+                    }
+                }
             }
         }
     }
@@ -1465,7 +1325,10 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
         if
         (
             is_infinite(c1) || is_infinite(c2)
-         || (!c1->hasInternalPoint() && !c2->hasInternalPoint())
+         || (
+                !c1->internalOrBoundaryDualVertex()
+             || !c2->internalOrBoundaryDualVertex()
+            )
          || !c1->real() || !c2->real()
         )
         {
@@ -1475,19 +1338,24 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
 //        Foam::point endPt = 0.5*(c1->dual() + c2->dual());
         Foam::point endPt = c1->dual();
 
-        if
-        (
-            magSqr(vert - c1->dual())
-          < magSqr(vert - c2->dual())
-        )
+        if (magSqr(vert - c1->dual()) < magSqr(vert - c2->dual()))
         {
             endPt = c2->dual();
         }
 
+        if
+        (
+            magSqr(vert - endPt)
+          > magSqr(geometryToConformTo().globalBounds().mag())
+        )
+        {
+            continue;
+        }
+
         pointIndexHit surfHit;
         label hitSurface;
 
-        geometryToConformTo_.findSurfaceAnyIntersection
+        geometryToConformTo_.findSurfaceNearestIntersection
         (
             vert,
             endPt,
@@ -1507,19 +1375,46 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
 
             const vector& n = norm[0];
 
-            scalar normalIncursionDistance = (endPt - surfHit.hitPoint()) & n;
+            scalar normalIncursionDistance
+            (
+                (endPt - surfHit.hitPoint()) & n
+            );
 
             if (normalIncursionDistance < minIncursionDistance)
             {
-                surfHitLargest = surfHit;
-                hitSurfaceLargest = hitSurface;
+                const plane p(surfHit.hitPoint(), n);
 
-                minIncursionDistance = normalIncursionDistance;
+                const plane::ray r(endPt, n);
 
-                // Info<< nl << "# Incursion: " << endl;
-                // meshTools::writeOBJ(Info, vert);
-                // meshTools::writeOBJ(Info, edgeMid);
-                // Info<< "l Na Nb" << endl;
+                const scalar d = p.normalIntersect(r);
+
+                Foam::point newPoint = endPt + d*n;
+
+                pointIndexHitAndFeature info;
+                geometryToConformTo_.findSurfaceNearest
+                (
+                    newPoint,
+                    4.0*magSqr(newPoint - endPt),
+                    info.first(),
+                    info.second()
+                );
+
+                if (info.first().hit())
+                {
+                    if
+                    (
+                        surfaceLocationConformsToInside
+                        (
+                            pointIndexHitAndFeature(info.first(), info.second())
+                        )
+                    )
+                    {
+                        surfHitLargest = info.first();
+                        hitSurfaceLargest = info.second();
+
+                        minIncursionDistance = normalIncursionDistance;
+                    }
+                }
             }
         }
     }
@@ -2219,6 +2114,9 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                             // feature edge, give control to edge control points
                             // instead, this will prevent "pits" forming.
 
+                            // Allow if different surfaces
+
+
                             keepSurfacePoint = false;
                         }
 
@@ -2298,11 +2196,13 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
             surfaceHits.append(surfHitI);
             appendToSurfacePtTree(surfPt);
             surfaceToTreeShape.append(existingSurfacePtLocations_.size() - 1);
+
+//            addedPoints.write(surfPt);
+        }
+        else
+        {
+//            removedPoints.write(surfPt);
         }
-//        else
-//        {
-//            surfaceToTreeShape.remove();
-//        }
     }
 }
 
@@ -2338,7 +2238,9 @@ void Foam::conformalVoronoiMesh::storeSurfaceConformation()
                 Vb
                 (
                     vit->point(),
-                    vit->type()
+                    vit->index(),
+                    vit->type(),
+                    Pstream::myProcNo()
                 )
             );
         }
@@ -2362,24 +2264,40 @@ void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
 {
     Info<< nl << "Reinserting stored surface conformation" << endl;
 
-    const label preReinsertionSize(number_of_vertices());
+    Map<label> oldToNewIndices =
+        insertPointPairs(surfaceConformationVertices_, true, true);
 
-    // It is assumed that the stored surface conformation is on the correct
-    // processor and does not need distributed
-    rangeInsertWithInfo
+    ptPairs_.reIndex(oldToNewIndices);
+
+    PackedBoolList selectedElems(surfaceConformationVertices_.size(), true);
+
+    forAll(surfaceConformationVertices_, vI)
+    {
+        Vb& v = surfaceConformationVertices_[vI];
+        label& vIndex = v.index();
+
+        Map<label>::const_iterator iter = oldToNewIndices.find(vIndex);
+
+        if (iter != oldToNewIndices.end())
+        {
+            const label newIndex = iter();
+
+            if (newIndex != -1)
+            {
+                vIndex = newIndex;
+            }
+            else
+            {
+                selectedElems[vI] = false;
+            }
+        }
+    }
+
+    inplaceSubset<PackedBoolList, List<Vb> >
     (
-        surfaceConformationVertices_.begin(),
-        surfaceConformationVertices_.end(),
-        true
+        selectedElems,
+        surfaceConformationVertices_
     );
-
-    const label nInserted = label(number_of_vertices()) - preReinsertionSize;
-    const label nFailed = surfaceConformationVertices_.size() - nInserted;
-
-    Info<< "    " << returnReduce(nInserted, sumOp<label>())
-        << " points reinserted, failed to insert "
-        << returnReduce(nFailed, sumOp<label>())
-        << endl;
 }
 
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
index 980064fb48..f0c4c7d746 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
@@ -29,6 +29,7 @@ License
 #include "tetrahedron.H"
 #include "const_circulator.H"
 #include "DelaunayMeshTools.H"
+#include "OBJstream.H"
 
 using namespace Foam::vectorTools;
 
@@ -225,21 +226,48 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
         vector masterPtVec(normalDir + nextNormalDir);
         masterPtVec /= mag(masterPtVec) + SMALL;
 
-        if (((normalDir ^ nextNormalDir) & edDir) < SMALL)
+        if
+        (
+            ((normalDir ^ nextNormalDir) & edDir) < SMALL
+         || mag(masterPtVec) < SMALL
+        )
         {
 //            Info<< "    IGNORE REGION" << endl;
             addedMasterPreviously = false;
 
-            if (circ.size() == 2 && mag((normal & nextNormal) - 1) < SMALL)
+            if
+            (
+                circ.size() == 2
+             && mag((normal & nextNormal) - 1) < SMALL
+            )
             {
-                // Add an extra point
-
                 const vector n = 0.5*(normal + nextNormal);
 
-                vector s = ppDist*(edDir ^ n);
+                const vector s = ppDist*(edDir ^ n);
 
-                createPointPair(ppDist, edgePt + s, n, pts);
-                createPointPair(ppDist, edgePt - s, n, pts);
+                if (volType == extendedFeatureEdgeMesh::BOTH)
+                {
+                    createBafflePointPair(ppDist, edgePt + s, n, true, pts);
+                    createBafflePointPair(ppDist, edgePt - s, n, true, pts);
+                }
+                else
+                {
+                    WarningIn
+                    (
+                        "Foam::conformalVoronoiMesh::"
+                        "createEdgePointGroupByCirculating"
+                        "("
+                        "    const extendedFeatureEdgeMesh&,"
+                        "    const pointIndexHit&,"
+                        "    DynamicList<Vb>&"
+                        ")"
+                    )   << "Failed to insert flat/open edge as volType is "
+                        << extendedFeatureEdgeMesh::sideVolumeTypeNames_
+                           [
+                               volType
+                           ]
+                        << endl;
+                }
 
                 break;
             }
@@ -247,34 +275,6 @@ void Foam::conformalVoronoiMesh::createEdgePointGroupByCirculating
             continue;
         }
 
-        if (mag(masterPtVec) < SMALL)
-        {
-            if (circ.size() == 2)
-            {
-                // Add an extra point
-
-                // Average normal to remove any bias to one side, although as
-                // it is a flat edge, the normals should be essentially the same
-                const vector n = 0.5*(normal + nextNormal);
-
-                // Direction along the surface to the control point, sense of
-                // direction not important, as +s and -s can be used because it
-                // is a flat edge
-                vector s = ppDist*(edDir ^ n);
-
-                createPointPair(ppDist, edgePt + s, n, pts);
-                createPointPair(ppDist, edgePt - s, n, pts);
-
-                break;
-            }
-            else
-            {
-//                Info<< "    IGNORE REGION" << endl;
-                addedMasterPreviously = false;
-                continue;
-            }
-        }
-
         const Foam::point masterPt = edgePt + ppDist*masterPtVec;
 
         // Check that region is inside or outside
@@ -490,11 +490,19 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
 
     const vectorField& feNormals = feMesh.normals();
     const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+    const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
+        feMesh.normalVolumeTypes();
 
     // As this is an external edge, there are two normals by definition
     const vector& nA = feNormals[edNormalIs[0]];
     const vector& nB = feNormals[edNormalIs[1]];
 
+    const extendedFeatureEdgeMesh::sideVolumeType& volTypeA =
+        normalVolumeTypes[edNormalIs[0]];
+
+    const extendedFeatureEdgeMesh::sideVolumeType& volTypeB =
+        normalVolumeTypes[edNormalIs[1]];
+
     if (areParallel(nA, nB))
     {
         // The normals are nearly parallel, so this is too sharp a feature to
@@ -524,13 +532,6 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
     // Convex. So refPt will be inside domain and hence a master point
     Foam::point refPt = edgePt - ppDist*refVec;
 
-    // Result when the points are eventually inserted.
-    // Add number_of_vertices() at insertion of first vertex to all numbers:
-    // pt           index type
-    // refPt        0     1
-    // reflectedA   1     0
-    // reflectedB   2     0
-
     // Insert the master point pairing the the first slave
 
     if (!geometryToConformTo_.inside(refPt))
@@ -559,7 +560,11 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
         (
             reflectedA,
             vertexCount() + pts.size(),
-            Vb::vtExternalFeatureEdge,
+            (
+                volTypeA == extendedFeatureEdgeMesh::BOTH
+              ? Vb::vtInternalFeatureEdge
+              : Vb::vtExternalFeatureEdge
+            ),
             Pstream::myProcNo()
         )
     );
@@ -571,10 +576,26 @@ void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
         (
             reflectedB,
             vertexCount() + pts.size(),
-            Vb::vtExternalFeatureEdge,
+            (
+                volTypeB == extendedFeatureEdgeMesh::BOTH
+              ? Vb::vtInternalFeatureEdge
+              : Vb::vtExternalFeatureEdge
+            ),
             Pstream::myProcNo()
         )
     );
+
+    ptPairs_.addPointPair
+    (
+        pts[pts.size() - 3].index(),
+        pts[pts.size() - 1].index()
+    );
+
+    ptPairs_.addPointPair
+    (
+        pts[pts.size() - 3].index(),
+        pts[pts.size() - 2].index()
+    );
 }
 
 
@@ -591,11 +612,19 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
 
     const vectorField& feNormals = feMesh.normals();
     const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+    const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
+        feMesh.normalVolumeTypes();
 
     // As this is an external edge, there are two normals by definition
     const vector& nA = feNormals[edNormalIs[0]];
     const vector& nB = feNormals[edNormalIs[1]];
 
+    const extendedFeatureEdgeMesh::sideVolumeType& volTypeA =
+        normalVolumeTypes[edNormalIs[0]];
+
+//    const extendedFeatureEdgeMesh::sideVolumeType& volTypeB =
+//        normalVolumeTypes[edNormalIs[1]];
+
     if (areParallel(nA, nB))
     {
         // The normals are nearly parallel, so this is too sharp a feature to
@@ -705,14 +734,30 @@ void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
         (
             reflMasterPt,
             vertexCount() + pts.size(),
-            Vb::vtExternalFeatureEdge,
+            (
+                volTypeA == extendedFeatureEdgeMesh::BOTH
+              ? Vb::vtInternalFeatureEdge
+              : Vb::vtExternalFeatureEdge
+            ),
             Pstream::myProcNo()
         )
     );
 
+    ptPairs_.addPointPair
+    (
+        pts[pts.size() - 2].index(),
+        pts[pts.size() - 1].index()
+    );
+
+    ptPairs_.addPointPair
+    (
+        pts[pts.size() - 3].index(),
+        pts[pts.size() - 1].index()
+    );
+
     if (nAddPoints == 1)
     {
-        // One additinal point is the reflection of the slave point,
+        // One additional point is the reflection of the slave point,
         // i.e. the original reference point
         pts.append
         (
@@ -767,6 +812,8 @@ void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
 
     const vectorField& feNormals = feMesh.normals();
     const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+    const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
+        feMesh.normalVolumeTypes();
 
     // As this is a flat edge, there are two normals by definition
     const vector& nA = feNormals[edNormalIs[0]];
@@ -781,8 +828,21 @@ void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
     // is a flat edge
     vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
 
-    createPointPair(ppDist, edgePt + s, n, pts);
-    createPointPair(ppDist, edgePt - s, n, pts);
+    if (normalVolumeTypes[edNormalIs[0]] == extendedFeatureEdgeMesh::OUTSIDE)
+    {
+        createPointPair(ppDist, edgePt + s, -n, true, pts);
+        createPointPair(ppDist, edgePt - s, -n, true, pts);
+    }
+    else if (normalVolumeTypes[edNormalIs[0]] == extendedFeatureEdgeMesh::BOTH)
+    {
+        createBafflePointPair(ppDist, edgePt + s, n, true, pts);
+        createBafflePointPair(ppDist, edgePt - s, n, true, pts);
+    }
+    else
+    {
+        createPointPair(ppDist, edgePt + s, n, true, pts);
+        createPointPair(ppDist, edgePt - s, n, true, pts);
+    }
 }
 
 
@@ -814,16 +874,23 @@ void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
 
         plane facePlane(edgePt, n);
 
-        Foam::point pt1 = edgePt + s + ppDist*n;
-        Foam::point pt2 = edgePt - s + ppDist*n;
+        const label initialPtsSize = pts.size();
 
-        Foam::point pt3 = facePlane.mirror(pt1);
-        Foam::point pt4 = facePlane.mirror(pt2);
+        if
+        (
+            !geometryToConformTo_.inside(edgePt)
+        )
+        {
+            return;
+        }
 
-        pts.append(Vb(pt1, Vb::vtInternalSurface));
-        pts.append(Vb(pt2, Vb::vtInternalSurface));
-        pts.append(Vb(pt3, Vb::vtInternalSurface));
-        pts.append(Vb(pt4, Vb::vtInternalSurface));
+        createBafflePointPair(ppDist, edgePt - s, n, true, pts);
+        createBafflePointPair(ppDist, edgePt + s, n, false, pts);
+
+        for (label ptI = initialPtsSize; ptI < pts.size(); ++ptI)
+        {
+            pts[ptI].type() = Vb::vtInternalFeatureEdge;
+        }
     }
     else
     {
@@ -846,26 +913,245 @@ void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
 
     const scalar ppDist = pointPairDistance(edgePt);
 
+    const vector edDir = feMesh.edgeDirections()[edHit.index()];
+
     const vectorField& feNormals = feMesh.normals();
     const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
+    const labelList& normalDirs = feMesh.normalDirections()[edHit.index()];
 
-    Info<< edNormalIs.size() << endl;
+    const List<extendedFeatureEdgeMesh::sideVolumeType>& normalVolumeTypes =
+        feMesh.normalVolumeTypes();
 
-    // As this is a flat edge, there are two normals by definition
-    const vector& nA = feNormals[edNormalIs[0]];
-    const vector& nB = feNormals[edNormalIs[1]];
+    labelList nNormalTypes(4, 0);
 
-    // Average normal to remove any bias to one side, although as this
-    // is a flat edge, the normals should be essentially the same
-    const vector n = 0.5*(nA + nB);
+    forAll(edNormalIs, edgeNormalI)
+    {
+        const extendedFeatureEdgeMesh::sideVolumeType sType =
+            normalVolumeTypes[edNormalIs[edgeNormalI]];
 
-    // Direction along the surface to the control point, sense of edge
-    // direction not important, as +s and -s can be used because this
-    // is a flat edge
-    vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
+        nNormalTypes[sType]++;
+    }
 
-    createPointPair(ppDist, edgePt + s, n, pts);
-    createPointPair(ppDist, edgePt - s, n, pts);
+    if (nNormalTypes[extendedFeatureEdgeMesh::BOTH] == 4)
+    {
+        label masterEdgeNormalIndex = -1;
+
+        forAll(edNormalIs, edgeNormalI)
+        {
+            const extendedFeatureEdgeMesh::sideVolumeType sType =
+                normalVolumeTypes[edNormalIs[edgeNormalI]];
+
+            if (sType == extendedFeatureEdgeMesh::BOTH)
+            {
+                masterEdgeNormalIndex = edgeNormalI;
+                break;
+            }
+        }
+
+        const vector& n = feNormals[edNormalIs[masterEdgeNormalIndex]];
+
+        label nDir = normalDirs[masterEdgeNormalIndex];
+
+        vector normalDir =
+            (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
+        normalDir *= nDir/mag(normalDir);
+
+        Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*n;
+        Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*n;
+
+        plane plane3(edgePt, normalDir);
+
+        Foam::point pt3 = plane3.mirror(pt1);
+        Foam::point pt4 = plane3.mirror(pt2);
+
+        pts.append
+        (
+            Vb
+            (
+                pt1,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+        pts.append
+        (
+            Vb
+            (
+                pt2,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 2].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                pt3,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 3].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                pt4,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 3].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 2].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+    }
+    else if
+    (
+        nNormalTypes[extendedFeatureEdgeMesh::BOTH] == 1
+     && nNormalTypes[extendedFeatureEdgeMesh::INSIDE] == 2
+    )
+    {
+        label masterEdgeNormalIndex = -1;
+
+        forAll(edNormalIs, edgeNormalI)
+        {
+            const extendedFeatureEdgeMesh::sideVolumeType sType =
+                normalVolumeTypes[edNormalIs[edgeNormalI]];
+
+            if (sType == extendedFeatureEdgeMesh::BOTH)
+            {
+                masterEdgeNormalIndex = edgeNormalI;
+                break;
+            }
+        }
+
+        const vector& n = feNormals[edNormalIs[masterEdgeNormalIndex]];
+
+        label nDir = normalDirs[masterEdgeNormalIndex];
+
+        vector normalDir =
+            (feNormals[edNormalIs[masterEdgeNormalIndex]] ^ edDir);
+        normalDir *= nDir/mag(normalDir);
+
+        const label nextNormalI =
+            (masterEdgeNormalIndex + 1) % edNormalIs.size();
+        if ((normalDir & feNormals[edNormalIs[nextNormalI]]) > 0)
+        {
+            normalDir *= -1;
+        }
+
+        Foam::point pt1 = edgePt + ppDist*normalDir + ppDist*n;
+        Foam::point pt2 = edgePt + ppDist*normalDir - ppDist*n;
+
+        plane plane3(edgePt, normalDir);
+
+        Foam::point pt3 = plane3.mirror(pt1);
+        Foam::point pt4 = plane3.mirror(pt2);
+
+        pts.append
+        (
+            Vb
+            (
+                pt1,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+        pts.append
+        (
+            Vb
+            (
+                pt2,
+                vertexCount() + pts.size(),
+                Vb::vtInternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 2].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                pt3,
+                vertexCount() + pts.size(),
+                Vb::vtExternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 3].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+
+        pts.append
+        (
+            Vb
+            (
+                pt4,
+                vertexCount() + pts.size(),
+                Vb::vtExternalSurface,
+                Pstream::myProcNo()
+            )
+        );
+
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 3].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+    }
+
+
+//    // As this is a flat edge, there are two normals by definition
+//    const vector& nA = feNormals[edNormalIs[0]];
+//    const vector& nB = feNormals[edNormalIs[1]];
+//
+//    // Average normal to remove any bias to one side, although as this
+//    // is a flat edge, the normals should be essentially the same
+//    const vector n = 0.5*(nA + nB);
+//
+//    // Direction along the surface to the control point, sense of edge
+//    // direction not important, as +s and -s can be used because this
+//    // is a flat edge
+//    vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
+//
+//    createBafflePointPair(ppDist, edgePt + s, n, true, pts);
+//    createBafflePointPair(ppDist, edgePt - s, n, true, pts);
 }
 
 
@@ -874,6 +1160,23 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints(bool distribute)
     Info<< nl
         << "Inserting feature points" << endl;
 
+
+//    Map<label> oldToNewIndices = insertPoints(featureVertices_, distribute);
+
+//    featurePtPairs_.reIndex(oldToNewIndices);
+
+//    forAll(featureVertices_, vI)
+//    {
+//        Map<label>::const_iterator iter =
+//            oldToNewIndices.find(featureVertices_[vI].index());
+
+//        if (iter != oldToNewIndices.end())
+//        {
+//            featureVertices_[vI].index() = iter();
+//        }
+//    }
+
+
     const label preFeaturePointSize(number_of_vertices());
 
     if (Pstream::parRun() && distribute)
@@ -884,7 +1187,7 @@ void Foam::conformalVoronoiMesh::insertFeaturePoints(bool distribute)
     const List<Vb>& ftPtVertices = ftPtConformer_.featurePointVertices();
 
     // Insert the created points directly as already distributed.
-    this->DelaunayMesh<Delaunay>::insertPoints(ftPtVertices);
+    this->DelaunayMesh<Delaunay>::insertPoints(ftPtVertices, true);
 
     label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
     reduce(nFeatureVertices, sumOp<label>());
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index e9e03a6330..7d35594258 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -225,6 +225,7 @@ inline void Foam::conformalVoronoiMesh::createPointPair
     const scalar ppDist,
     const Foam::point& surfPt,
     const vector& n,
+    const bool ptPair,
     DynamicList<Vb>& pts
 ) const
 {
@@ -260,14 +261,14 @@ inline void Foam::conformalVoronoiMesh::createPointPair
             )
         );
 
-//        if (ptPair)
-//        {
-//            ptPairs_.addPointPair
-//            (
-//                pts[pts.size() - 2].index(),
-//                pts[pts.size() - 1].index() // external 0 -> slave
-//            );
-//        }
+        if (ptPair)
+        {
+            ptPairs_.addPointPair
+            (
+                pts[pts.size() - 2].index(),
+                pts[pts.size() - 1].index() // external 0 -> slave
+            );
+        }
     }
 //    else
 //    {
@@ -305,6 +306,7 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
     const scalar ppDist,
     const Foam::point& surfPt,
     const vector& n,
+    const bool ptPair,
     DynamicList<Vb>& pts
 ) const
 {
@@ -315,8 +317,8 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
         Vb
         (
             surfPt - ppDistn,
-            vertexCount() + 1 + pts.size(),
-            Vb::vtInternalSurface,
+            vertexCount() + pts.size(),
+            Vb::vtInternalSurfaceBaffle,
             Pstream::myProcNo()
         )
     );
@@ -326,11 +328,20 @@ inline void Foam::conformalVoronoiMesh::createBafflePointPair
         Vb
         (
             surfPt + ppDistn,
-            vertexCount() + 1 + pts.size(),
-            Vb::vtInternalSurface,
+            vertexCount() + pts.size(),
+            Vb::vtExternalSurfaceBaffle,
             Pstream::myProcNo()
         )
     );
+
+    if (ptPair)
+    {
+        ptPairs_.addPointPair
+        (
+            pts[pts.size() - 2].index(), // external 0 -> slave
+            pts[pts.size() - 1].index()
+        );
+    }
 }
 
 
@@ -341,8 +352,8 @@ inline bool Foam::conformalVoronoiMesh::internalPointIsInside
 {
     if
     (
-        !geometryToConformTo_.inside(pt)
-     || !geometryToConformTo_.globalBounds().contains(pt)
+        !geometryToConformTo_.globalBounds().contains(pt)
+     || !geometryToConformTo_.inside(pt)
     )
     {
         return false;
@@ -352,76 +363,6 @@ inline bool Foam::conformalVoronoiMesh::internalPointIsInside
 }
 
 
-inline bool Foam::conformalVoronoiMesh::isPointPair
-(
-    const Vertex_handle& vA,
-    const Vertex_handle& vB
-) const
-{
-    // Want to do this topologically, but problem if vertices are redistributed
-    // so that one of the point pair is one processor and the other is on
-    // another.
-
-    const Foam::point& ptA = topoint(vA->point());
-    const Foam::point& ptB = topoint(vB->point());
-
-    if
-    (
-        (
-            vA->type() == Vb::vtInternalSurface
-         && vB->type() == Vb::vtExternalSurface
-        )
-     ||
-        (
-            vB->type() == Vb::vtInternalSurface
-         && vA->type() == Vb::vtExternalSurface
-        )
-     ||
-        (
-            vA->type() == Vb::vtInternalFeatureEdge
-         && vB->type() == Vb::vtExternalFeatureEdge
-        )
-     ||
-        (
-            vB->type() == Vb::vtInternalFeatureEdge
-         && vA->type() == Vb::vtExternalFeatureEdge
-        )
-     ||
-        (
-            vA->type() == Vb::vtInternalSurface
-         && vB->type() == Vb::vtExternalFeatureEdge
-        )
-     ||
-        (
-            vB->type() == Vb::vtInternalSurface
-         && vA->type() == Vb::vtExternalFeatureEdge
-        )
-     ||
-        (
-            vA->type() == Vb::vtExternalSurface
-         && vB->type() == Vb::vtInternalFeatureEdge
-        )
-     ||
-        (
-            vB->type() == Vb::vtExternalSurface
-         && vA->type() == Vb::vtInternalFeatureEdge
-        )
-    )
-    {
-        const scalar distSqr = magSqr(ptA - ptB);
-
-        const scalar ppDistSqr = sqr(2*pointPairDistance(0.5*(ptA + ptB)));
-
-        if (distSqr > 1.001*ppDistSqr)
-        {
-            return false;
-        }
-    }
-
-    return true;
-}
-
-
 inline bool Foam::conformalVoronoiMesh::isBoundaryDualFace
 (
     const Delaunay::Finite_edges_iterator& eit
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index d733b33e35..b0da21cbbb 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -36,6 +36,8 @@ License
 #include "indexedVertexOps.H"
 #include "DelaunayMeshTools.H"
 #include "syncTools.H"
+#include "faceSet.H"
+#include "OBJstream.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -114,7 +116,6 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
         wordList patchNames;
         PtrList<dictionary> patchDicts;
         pointField cellCentres;
-
         PackedBoolList boundaryFacesToRemove;
 
         calcDualMesh
@@ -789,7 +790,7 @@ void Foam::conformalVoronoiMesh::writeMesh
     const wordList& patchNames,
     const PtrList<dictionary>& patchDicts,
     const pointField& cellCentres,
-    const PackedBoolList& boundaryFacesToRemove
+    PackedBoolList& boundaryFacesToRemove
 ) const
 {
     if (foamyHexMeshControls().objOutput())
@@ -912,42 +913,84 @@ void Foam::conformalVoronoiMesh::writeMesh
 
 
 
-    // Add indirectPatchFaces to a face zone
+    Info<< indent << "Add pointZones" << endl;
     {
-        labelList addr(boundaryFacesToRemove.count());
-        label count = 0;
+        label sz = mesh.pointZones().size();
 
-        forAll(boundaryFacesToRemove, faceI)
+        DynamicList<label> bPts(boundaryPts.size());
+
+        forAll(dualMeshPointTypeNames_, typeI)
         {
-            if
-            (
-                boundaryFacesToRemove[faceI]
-             && mesh.faceZones().whichZone(faceI) == -1
-            )
+            forAll(boundaryPts, ptI)
             {
-                addr[count++] = faceI;
+                const label& bPtType = boundaryPts[ptI];
+
+                if (bPtType == typeI)
+                {
+                    bPts.append(ptI);
+                }
             }
-        }
 
-        addr.setSize(count);
+//            syncTools::syncPointList(mesh, bPts, maxEqOp<label>(), -1);
 
-        label sz = mesh.faceZones().size();
-        boolList flip(addr.size(), false);
-        mesh.faceZones().setSize(sz + 1);
-        mesh.faceZones().set
-        (
-            sz,
-            new faceZone
+            Info<< incrIndent << indent
+                << "Adding " << bPts.size()
+                << " points of type " << dualMeshPointTypeNames_.words()[typeI]
+                << decrIndent << endl;
+
+            mesh.pointZones().append
             (
-                "indirectPatchFaces",
-                addr,
-                flip,
-                sz,
-                mesh.faceZones()
-            )
-        );
+                new pointZone
+                (
+                    dualMeshPointTypeNames_.words()[typeI],
+                    bPts,
+                    sz + typeI,
+                    mesh.pointZones()
+                )
+            );
+
+            bPts.clear();
+        }
     }
 
+
+
+    // Add indirectPatchFaces to a face zone
+    Info<< indent << "Adding indirect patch faces set" << endl;
+
+    syncTools::syncFaceList
+    (
+        mesh,
+        boundaryFacesToRemove,
+        orOp<bool>()
+    );
+
+    labelList addr(boundaryFacesToRemove.count());
+    label count = 0;
+
+    forAll(boundaryFacesToRemove, faceI)
+    {
+        if (boundaryFacesToRemove[faceI])
+        {
+            addr[count++] = faceI;
+        }
+    }
+
+    addr.setSize(count);
+
+    faceSet indirectPatchFaces
+    (
+        mesh,
+        "indirectPatchFaces",
+        addr,
+        IOobject::AUTO_WRITE
+    );
+
+    indirectPatchFaces.sync(mesh);
+
+
+    Info<< decrIndent;
+
     timeCheck("Before fvMesh filtering");
 
     autoPtr<polyMeshFilter> meshFilter;
@@ -966,9 +1009,11 @@ void Foam::conformalVoronoiMesh::writeMesh
         {
             const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
 
-            polyTopoChange meshMod(newMesh);
+            polyTopoChange meshMod(newMesh());
 
-            meshMod.changeMesh(mesh, false);
+            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+            polyMeshFilter::copySets(newMesh(), mesh);
         }
     }
 
@@ -1001,9 +1046,11 @@ void Foam::conformalVoronoiMesh::writeMesh
         {
             const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
 
-            polyTopoChange meshMod(newMesh);
+            polyTopoChange meshMod(newMesh());
 
-            meshMod.changeMesh(mesh, false);
+            autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+            polyMeshFilter::copySets(newMesh(), mesh);
         }
     }
 
@@ -1068,7 +1115,7 @@ void Foam::conformalVoronoiMesh::writeMesh
 
 //    writeCellCentres(mesh);
 
-//    findRemainingProtrusionSet(mesh);
+    findRemainingProtrusionSet(mesh);
 }
 
 
@@ -1375,4 +1422,33 @@ Foam::labelHashSet Foam::conformalVoronoiMesh::findRemainingProtrusionSet
 }
 
 
+void Foam::conformalVoronoiMesh::writePointPairs
+(
+    const fileName& fName
+) const
+{
+    OBJstream os(fName);
+
+    for
+    (
+        Delaunay::Finite_edges_iterator eit = finite_edges_begin();
+        eit != finite_edges_end();
+        ++eit
+    )
+    {
+        Cell_handle c = eit->first;
+        Vertex_handle vA = c->vertex(eit->second);
+        Vertex_handle vB = c->vertex(eit->third);
+
+        if (ptPairs_.isPointPair(vA, vB))
+        {
+            os.write
+            (
+                linePointRef(topoint(vA->point()), topoint(vB->point()))
+            );
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.C
index 542f84afe5..6501d72758 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.C
@@ -146,7 +146,7 @@ void Foam::featurePointConformer::addMasterAndSlavePoints
             )
         );
 
-//        const label masterIndex = pts[pts.size() - 1].index();
+        const label masterIndex = pts.last().index();
 
         //meshTools::writeOBJ(strMasters, masterPt);
 
@@ -180,6 +180,11 @@ void Foam::featurePointConformer::addMasterAndSlavePoints
                 )
             );
 
+            ftPtPairs_.addPointPair
+            (
+                masterIndex,
+                pts.last().index()
+            );
 
             //meshTools::writeOBJ(strSlaves, slavePt);
         }
@@ -525,7 +530,8 @@ Foam::featurePointConformer::featurePointConformer
     foamyHexMesh_(foamyHexMesh),
     foamyHexMeshControls_(foamyHexMesh.foamyHexMeshControls()),
     geometryToConformTo_(foamyHexMesh.geometryToConformTo()),
-    featurePointVertices_()
+    featurePointVertices_(),
+    ftPtPairs_(foamyHexMesh)
 {
     Info<< nl
         << "Conforming to feature points" << endl;
@@ -598,6 +604,8 @@ void Foam::featurePointConformer::distribute
     {
         featurePointVertices_[vI].procIndex() = Pstream::myProcNo();
     }
+
+    // Also update the feature point pairs
 }
 
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.H
index d4d52c1fff..ddda0959e8 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformer.H
@@ -44,6 +44,7 @@ SourceFiles
 #include "DynamicList.H"
 #include "List.H"
 #include "extendedFeatureEdgeMesh.H"
+#include "pointPairs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -83,6 +84,9 @@ class featurePointConformer
         //  triangulation clear.
         List<Vb> featurePointVertices_;
 
+        //-
+        mutable pointPairs<Delaunay> ftPtPairs_;
+
 
     // Private Member Functions
 
@@ -163,6 +167,9 @@ public:
             //  triangulation.
             inline const List<Vb>& featurePointVertices() const;
 
+            //- Return the feature point pair table
+            inline const pointPairs<Delaunay>& featurePointPairs() const;
+
 
         // Edit
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerI.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerI.H
index 406ed3f000..0666f274ba 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerI.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerI.H
@@ -28,5 +28,11 @@ const Foam::List<Vb>& Foam::featurePointConformer::featurePointVertices() const
     return featurePointVertices_;
 }
 
+const Foam::pointPairs<Delaunay>&
+Foam::featurePointConformer::featurePointPairs() const
+{
+    return ftPtPairs_;
+}
+
 
 // ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerSpecialisations.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerSpecialisations.C
index 5ddccc0d68..b12f56b909 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerSpecialisations.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/featurePointConformer/featurePointConformerSpecialisations.C
@@ -168,9 +168,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(internalPtA, Vb::vtInternalFeaturePoint)
+            Vb
+            (
+                internalPtA,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtInternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        const label internalPtAIndex(pts.last().index());
+
         const Foam::point internalPtB =
             concaveEdgeExternalPt
           - 2.0*planeB.distance(concaveEdgeExternalPt)
@@ -178,9 +186,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(internalPtB, Vb::vtInternalFeaturePoint)
+            Vb
+            (
+                internalPtB,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtInternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        const label internalPtBIndex(pts.last().index());
+
         // Add the external points
 
         Foam::point externalPtD;
@@ -288,7 +304,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
                         pts.append
                         (
-                            Vb(externalPtD, Vb::vtExternalFeaturePoint)
+                            Vb
+                            (
+                                externalPtD,
+                                foamyHexMesh_.vertexCount() + pts.size(),
+                                Vb::vtExternalFeaturePoint,
+                                Pstream::myProcNo()
+                            )
+                        );
+
+                        ftPtPairs_.addPointPair
+                        (
+                            internalPtAIndex,
+                            pts.last().index()
                         );
                     }
                 }
@@ -319,7 +347,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
                         pts.append
                         (
-                            Vb(externalPtE, Vb::vtExternalFeaturePoint)
+                            Vb
+                            (
+                                externalPtE,
+                                foamyHexMesh_.vertexCount() + pts.size(),
+                                Vb::vtExternalFeaturePoint,
+                                Pstream::myProcNo()
+                            )
+                        );
+
+                        ftPtPairs_.addPointPair
+                        (
+                            internalPtBIndex,
+                            pts.last().index()
                         );
                     }
                 }
@@ -328,9 +368,29 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(concaveEdgeExternalPt, Vb::vtExternalFeaturePoint)
+            Vb
+            (
+                concaveEdgeExternalPt,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtExternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        ftPtPairs_.addPointPair
+        (
+            internalPtBIndex,
+            pts.last().index()
+        );
+
+        ftPtPairs_.addPointPair
+        (
+            internalPtAIndex,
+            pts.last().index()
+        );
+
+        const label concaveEdgeExternalPtIndex(pts.last().index());
+
         const scalar totalAngle = radToDeg
         (
             constant::mathematical::pi
@@ -377,7 +437,21 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
             pts.append
             (
-                Vb(internalPtF, Vb::vtInternalFeaturePoint)
+                Vb
+                (
+                    internalPtF,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    Vb::vtInternalFeaturePoint,
+                    Pstream::myProcNo()
+                )
+            );
+
+            const label internalPtFIndex(pts.last().index());
+
+            ftPtPairs_.addPointPair
+            (
+                concaveEdgeExternalPtIndex,
+                pts.last().index()
             );
 
             const Foam::point externalPtG =
@@ -386,7 +460,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
             pts.append
             (
-                Vb(externalPtG, Vb::vtExternalFeaturePoint)
+                Vb
+                (
+                    externalPtG,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    Vb::vtExternalFeaturePoint,
+                    Pstream::myProcNo()
+                )
+            );
+
+            ftPtPairs_.addPointPair
+            (
+                internalPtFIndex,
+                pts.last().index()
             );
         }
 
@@ -520,9 +606,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(internalPtA, Vb::vtExternalFeaturePoint)
+            Vb
+            (
+                internalPtA,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtExternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        const label internalPtAIndex(pts.last().index());
+
         const Foam::point internalPtB =
             convexEdgeExternalPt
           + 2.0*planeB.distance(convexEdgeExternalPt)
@@ -530,9 +624,17 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(internalPtB, Vb::vtExternalFeaturePoint)
+            Vb
+            (
+                internalPtB,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtExternalFeaturePoint,
+                Pstream::myProcNo()
+            )
         );
 
+        const label internalPtBIndex(pts.last().index());
+
         // Add the internal points
 
         Foam::point externalPtD;
@@ -643,7 +745,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
                         pts.append
                         (
-                            Vb(externalPtD, Vb::vtInternalFeaturePoint)
+                            Vb
+                            (
+                                externalPtD,
+                                foamyHexMesh_.vertexCount() + pts.size(),
+                                Vb::vtInternalFeaturePoint,
+                                Pstream::myProcNo()
+                            )
+                        );
+
+                        ftPtPairs_.addPointPair
+                        (
+                            internalPtAIndex,
+                            pts.last().index()
                         );
                     }
                 }
@@ -674,7 +788,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
                         pts.append
                         (
-                            Vb(externalPtE, Vb::vtInternalFeaturePoint)
+                            Vb
+                            (
+                                externalPtE,
+                                foamyHexMesh_.vertexCount() + pts.size(),
+                                Vb::vtInternalFeaturePoint,
+                                Pstream::myProcNo()
+                            )
+                        );
+
+                        ftPtPairs_.addPointPair
+                        (
+                            internalPtBIndex,
+                            pts.last().index()
                         );
                     }
                 }
@@ -683,7 +809,25 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
         pts.append
         (
-            Vb(convexEdgeExternalPt, Vb::vtInternalFeaturePoint)
+            Vb
+            (
+                convexEdgeExternalPt,
+                foamyHexMesh_.vertexCount() + pts.size(),
+                Vb::vtInternalFeaturePoint,
+                Pstream::myProcNo()
+            )
+        );
+
+        ftPtPairs_.addPointPair
+        (
+            internalPtBIndex,
+            pts.last().index()
+        );
+
+        ftPtPairs_.addPointPair
+        (
+            internalPtAIndex,
+            pts.last().index()
         );
 
         const scalar totalAngle = radToDeg
@@ -732,7 +876,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
             pts.append
             (
-                Vb(internalPtF, Vb::vtExternalFeaturePoint)
+                Vb
+                (
+                    internalPtF,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    Vb::vtExternalFeaturePoint,
+                    Pstream::myProcNo()
+                )
+            );
+
+            ftPtPairs_.addPointPair
+            (
+                pts[pts.size() - 2].index(),
+                pts.last().index()
             );
 
             const Foam::point externalPtG =
@@ -741,7 +897,19 @@ bool Foam::featurePointConformer::createSpecialisedFeaturePoint
 
             pts.append
             (
-                Vb(externalPtG, Vb::vtInternalFeaturePoint)
+                Vb
+                (
+                    externalPtG,
+                    foamyHexMesh_.vertexCount() + pts.size(),
+                    Vb::vtInternalFeaturePoint,
+                    Pstream::myProcNo()
+                )
+            );
+
+            ftPtPairs_.addPointPair
+            (
+                pts[pts.size() - 2].index(),
+                pts.last().index()
             );
         }
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
index 15b155d85b..977f6fd7ff 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCell.H
@@ -193,6 +193,8 @@ public:
         //- Does the Dual vertex form part of a processor patch
         inline bool parallelDualVertex() const;
 
+        inline Foam::label vertexLowestProc() const;
+
         //- Using the globalIndex object, return a list of four (sorted) global
         //  Delaunay vertex indices that uniquely identify this tet in parallel
         inline Foam::tetCell vertexGlobalIndices
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
index 6c55a8eb30..f4890b50fb 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H
@@ -291,6 +291,23 @@ inline bool CGAL::indexedCell<Gt, Cb>::parallelDualVertex() const
 }
 
 
+template<class Gt, class Cb>
+inline Foam::label CGAL::indexedCell<Gt, Cb>::vertexLowestProc() const
+{
+    Foam::label lowestProc = -1;
+
+    for (int i = 0; i < 4; ++i)
+    {
+        if (this->vertex(i)->referred())
+        {
+            lowestProc = min(lowestProc, this->vertex(i)->procIndex());
+        }
+    }
+
+    return lowestProc;
+}
+
+
 template<class Gt, class Cb>
 inline Foam::tetCell CGAL::indexedCell<Gt, Cb>::vertexGlobalIndices
 (
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 43d9c77013..ed944d86ce 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -83,6 +83,13 @@ void Foam::conformationSurfaces::hasBoundedVolume
 
             const pointField& surfPts = triSurf.points();
 
+            Info<< "    Checking " << surface.name() << endl;
+
+            label nBaffles = 0;
+
+            Info<< "        Index = " << surfaces_[s] << endl;
+            Info<< "        Offset = " << regionOffset_[s] << endl;
+
             forAll(triSurf, sI)
             {
                 const label patchID =
@@ -98,7 +105,13 @@ void Foam::conformationSurfaces::hasBoundedVolume
                 {
                     sum += triSurf[sI].normal(surfPts);
                 }
+                else
+                {
+                    nBaffles++;
+                }
             }
+            Info<< "        has " << nBaffles << " baffles out of "
+                << triSurf.size() << " triangles" << endl;
 
             totalTriangles += triSurf.size();
         }
@@ -153,8 +166,9 @@ void Foam::conformationSurfaces::readFeatures
     {
         const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
 
-        Info<< "    features: " << surface.name() << " of type "
-            << surface.type() << endl;
+        Info<< "    features: " << surface.name()
+            << " of type " << surface.type()
+            << ", id: " << featureIndex << endl;
 
         autoPtr<searchableSurfaceFeatures> ssFeatures
         (
@@ -210,7 +224,8 @@ void Foam::conformationSurfaces::readFeatures
     {
         fileName feMeshName(featureDict.lookup("extendedFeatureEdgeMesh"));
 
-        Info<< "    features: " << feMeshName << endl;
+        Info<< "    features: " << feMeshName << ", id: " << featureIndex
+            << endl;
 
         features_.set
         (
@@ -716,7 +731,11 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
 
             const searchableSurface& surface(allGeometry_[surfaces_[s]]);
 
-            if (!surface.hasVolumeType())
+            if
+            (
+                !surface.hasVolumeType()
+             //&& surfaceVolumeTests[s][i] == volumeType::UNKNOWN
+            )
             {
                 pointField sample(1, samplePts[i]);
                 scalarField nearestDistSqr(1, GREAT);
@@ -733,7 +752,7 @@ Foam::Field<bool> Foam::conformationSurfaces::wellInOutSide
                 findSurfaceNearestIntersection
                 (
                     samplePts[i],
-                    info[0].rawPoint(),
+                    info[0].rawPoint() - 1e-3*mag(hitDir)*hitDir,
                     surfHit,
                     hitSurface
                 );
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
index 5649e61480..3552323d68 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
@@ -187,7 +187,7 @@ Foam::label Foam::autoDensity::recurseAndFill
     word recursionName
 ) const
 {
-    label maxDepth = 0;
+    label treeDepth = 0;
 
     for (direction i = 0; i < 8; i++)
     {
@@ -201,10 +201,10 @@ Foam::label Foam::autoDensity::recurseAndFill
         {
             if (levelLimit > 0)
             {
-                maxDepth =
+                treeDepth =
                     max
                     (
-                        maxDepth,
+                        treeDepth,
                         recurseAndFill
                         (
                             initialPoints,
@@ -229,10 +229,10 @@ Foam::label Foam::autoDensity::recurseAndFill
 
                 if (!fillBox(initialPoints, subBB, true))
                 {
-                    maxDepth =
+                    treeDepth =
                         max
                         (
-                            maxDepth,
+                            treeDepth,
                             recurseAndFill
                             (
                                 initialPoints,
@@ -259,10 +259,10 @@ Foam::label Foam::autoDensity::recurseAndFill
 
             if (!fillBox(initialPoints, subBB, false))
             {
-                maxDepth =
+                treeDepth =
                     max
                     (
-                        maxDepth,
+                        treeDepth,
                         recurseAndFill
                         (
                             initialPoints,
@@ -286,7 +286,7 @@ Foam::label Foam::autoDensity::recurseAndFill
         }
     }
 
-    return maxDepth + 1;
+    return treeDepth + 1;
 }
 
 
@@ -941,7 +941,7 @@ List<Vb::Point> autoDensity::initialPoints() const
         Pout<< "    Filling box " << hierBB << endl;
     }
 
-    label maxDepth = recurseAndFill
+    label treeDepth = recurseAndFill
     (
         initialPoints,
         hierBB,
@@ -966,7 +966,7 @@ List<Vb::Point> autoDensity::initialPoints() const
         << scalar(nInitialPoints)/scalar(max(globalTrialPoints_, 1))
         << " success rate" << nl
         << indent
-        << returnReduce(maxDepth, maxOp<label>())
+        << returnReduce(treeDepth, maxOp<label>())
         << " levels of recursion (maximum)"
         << decrIndent << decrIndent
         << endl;
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.C
new file mode 100644
index 0000000000..901d81ab33
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.C
@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "pointPairs.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class Triangulation>
+inline Foam::Pair<Foam::labelPair>
+Foam::pointPairs<Triangulation>::orderPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+) const
+{
+    return
+    (
+        (vA < vB)
+      ? Pair<labelPair>(vA, vB)
+      : Pair<labelPair>(vB, vA)
+    );
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::findPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+) const
+{
+    if (vA == vB)
+    {
+        return false;
+    }
+    else if (find(orderPointPair(vA, vB)) == end())
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::insertPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+)
+{
+    if (vA == vB)
+    {
+        return false;
+    }
+
+    return insert(orderPointPair(vA, vB));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Triangulation>
+Foam::pointPairs<Triangulation>::pointPairs(const Triangulation& triangulation)
+:
+    ptPairTable(),
+    triangulation_(triangulation)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Triangulation>
+Foam::pointPairs<Triangulation>::~pointPairs()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::addPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+)
+{
+    return insertPointPair(vA, vB);
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::addPointPair
+(
+    const labelPair& master,
+    const DynamicList<labelPair>& slaves
+)
+{
+    forAll(slaves, sI)
+    {
+        const labelPair& slave = slaves[sI];
+
+        addPointPair(master, slave);
+    }
+
+    return true;
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::addPointPair
+(
+    const label vA,
+    const label vB
+)
+{
+    const label procNo = Pstream::myProcNo();
+
+    labelPair a(vA, procNo);
+    labelPair b(vB, procNo);
+
+    return addPointPair(a, b);
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::isPointPair
+(
+    const Vertex_handle& vA,
+    const Vertex_handle& vB
+) const
+{
+    if (vA->boundaryPoint() && vB->boundaryPoint())
+    {
+        labelPair a(vA->index(), vA->procIndex());
+        labelPair b(vB->index(), vB->procIndex());
+
+        return findPointPair(a, b);
+    }
+
+    return false;
+}
+
+
+template<class Triangulation>
+inline bool Foam::pointPairs<Triangulation>::isPointPair
+(
+    const labelPair& vA,
+    const labelPair& vB
+) const
+{
+    return findPointPair(vA, vB);
+}
+
+
+template<class Triangulation>
+void Foam::pointPairs<Triangulation>::reIndex(const Map<label>& oldToNewIndices)
+{
+    pointPairs<Triangulation> newTable(triangulation_);
+
+    forAllConstIter(pointPairs, *this, iter)
+    {
+        Pair<labelPair> e = iter.key();
+
+        labelPair& start = e.first();
+        labelPair& end = e.second();
+
+        bool insert = true;
+
+        if (start.second() == Pstream::myProcNo())
+        {
+            Map<label>::const_iterator iter2 =
+                oldToNewIndices.find(start.first());
+
+            if (iter2 != oldToNewIndices.end())
+            {
+                if (iter2() != -1)
+                {
+                    start.first() = iter2();
+                }
+                else
+                {
+                    insert = false;
+                }
+            }
+        }
+
+        if (end.second() == Pstream::myProcNo())
+        {
+            Map<label>::const_iterator iter2 =
+                oldToNewIndices.find(end.first());
+
+            if (iter2 != oldToNewIndices.end())
+            {
+                if (iter2() != -1)
+                {
+                    end.first() = iter2();
+                }
+                else
+                {
+                    insert = false;
+                }
+            }
+        }
+
+        if (insert)
+        {
+            if (e.first() < e.second())
+            {
+                newTable.insert(e);
+            }
+            else if (e.first() > e.second())
+            {
+                newTable.insert(reverse(e));
+            }
+        }
+    }
+
+    this->transfer(newTable);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.H
new file mode 100644
index 0000000000..b5fcd413ee
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/pointPairs/pointPairs.H
@@ -0,0 +1,163 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    pointPairs
+
+Description
+    HashSet of unique edges. The edges are stored as a pair of pairs:
+
+        ( (local index, processor index) (local index, processor index) )
+
+    e.g.,
+
+        ( (0 1) (3 1) )
+        ( (0 2) (5 1) )
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointPairs_H
+#define pointPairs_H
+
+#include "labelPair.H"
+#include "HashSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+typedef HashSet
+<
+    Pair<labelPair>,
+    FixedList<labelPair, 2>::Hash<>
+> ptPairTable;
+
+/*---------------------------------------------------------------------------*\
+                         Class pointPairs Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Triangulation>
+class pointPairs
+:
+    public ptPairTable
+{
+    // Private typedefs
+
+        typedef typename Triangulation::Vertex_handle Vertex_handle;
+
+
+    // Private data
+
+        const Triangulation& triangulation_;
+
+
+    // Private Member Functions
+
+        inline Pair<labelPair> orderPointPair
+        (
+            const labelPair& vA,
+            const labelPair& vB
+        ) const;
+
+        inline bool insertPointPair
+        (
+            const labelPair& vA,
+            const labelPair& vB
+        );
+
+        inline bool findPointPair
+        (
+            const labelPair& vA,
+            const labelPair& vB
+        ) const;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from triangulation
+        pointPairs(const Triangulation& triangulation);
+
+
+    //- Destructor
+    ~pointPairs();
+
+
+    // Member Functions
+
+        // Access
+
+            inline bool isPointPair
+            (
+                const Vertex_handle& vA,
+                const Vertex_handle& vB
+            ) const;
+
+            inline bool isPointPair
+            (
+                const labelPair& vA,
+                const labelPair& vB
+            ) const;
+
+
+        // Edit
+
+            inline bool addPointPair
+            (
+                const labelPair& vA,
+                const labelPair& vB
+            );
+
+            inline bool addPointPair
+            (
+                const labelPair& master,
+                const DynamicList<labelPair>& slaves
+            );
+
+            inline bool addPointPair
+            (
+                const label vA,
+                const label vB
+            );
+
+            void reIndex(const Map<label>& oldToNewIndices);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "pointPairs.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C
index 78149fac16..76726ccc07 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchableBoxFeatures.C
@@ -104,11 +104,14 @@ Foam::searchableBoxFeatures::features() const
     edgeNormals[10][0] = 1; edgeNormals[10][1] = 3;
     edgeNormals[11][0] = 3; edgeNormals[11][1] = 0;
 
+    tmp<pointField> surfacePointsTmp(surface().points());
+    pointField& surfacePoints = surfacePointsTmp();
+
     forAll(edgeDirections, eI)
     {
         edgeDirections[eI] =
-            surface().points()()[treeBoundBox::edges[eI].end()]
-          - surface().points()()[treeBoundBox::edges[eI].start()];
+            surfacePoints[treeBoundBox::edges[eI].end()]
+          - surfacePoints[treeBoundBox::edges[eI].start()];
 
         normalDirections[eI] = labelList(2, 0);
         for (label j = 0; j < 2; ++j)
@@ -116,8 +119,8 @@ Foam::searchableBoxFeatures::features() const
             const vector cross =
                 (faceNormals[edgeNormals[eI][j]] ^ edgeDirections[eI]);
             const vector fC0tofE0 =
-                0.5*(max(surface().points()() + min(surface().points()())))
-              - surface().points()()[treeBoundBox::edges[eI].start()];
+                0.5*(max(surfacePoints + min(surfacePoints)))
+              - surfacePoints[treeBoundBox::edges[eI].start()];
 
             normalDirections[eI][j] =
                 (
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.C
new file mode 100644
index 0000000000..a61862c432
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.C
@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "searchablePlateFeatures.H"
+#include "addToRunTimeSelectionTable.H"
+#include "treeBoundBox.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(searchablePlateFeatures, 0);
+addToRunTimeSelectionTable
+(
+    searchableSurfaceFeatures,
+    searchablePlateFeatures,
+    dict
+);
+
+//! \cond - skip documentation : local scope only
+const Foam::label edgesArray[4][2] =
+{
+    {0, 1}, // 0
+    {0, 3},
+    {2, 1}, // 2
+    {2, 3}
+};
+//! \endcond
+
+const edgeList searchablePlateFeatures::edges(calcEdges(edgesArray));
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::edgeList Foam::searchablePlateFeatures::calcEdges
+(
+    const label edgesArray[4][2]
+)
+{
+    edgeList edges(4);
+    forAll(edges, edgeI)
+    {
+        edges[edgeI][0] = edgesArray[edgeI][0];
+        edges[edgeI][1] = edgesArray[edgeI][1];
+    }
+    return edges;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::searchablePlateFeatures::searchablePlateFeatures
+(
+    const searchableSurface& surface,
+    const dictionary& dict
+)
+:
+    searchableSurfaceFeatures(surface, dict),
+    mode_
+    (
+        extendedFeatureEdgeMesh::sideVolumeTypeNames_
+        [
+            dict.lookupOrDefault<word>("meshableSide", "inside")
+        ]
+    )
+{
+    Info<< indent
+        << "    Meshable region = "
+        << extendedFeatureEdgeMesh::sideVolumeTypeNames_[mode_]
+        << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::searchablePlateFeatures::~searchablePlateFeatures()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::extendedFeatureEdgeMesh>
+Foam::searchablePlateFeatures::features() const
+{
+    autoPtr<extendedFeatureEdgeMesh> features;
+
+    vectorField faceNormals(1);
+    surface().getNormal(List<pointIndexHit>(1), faceNormals);
+
+    vectorField edgeDirections(4);
+    labelListList normalDirections(4);
+
+    labelListList edgeNormals(4);
+    forAll(edgeNormals, eI)
+    {
+        edgeNormals[eI].setSize(2, 0);
+    }
+    edgeNormals[0][0] = 0; edgeNormals[0][1] = 0;
+    edgeNormals[1][0] = 0; edgeNormals[1][1] = 0;
+    edgeNormals[2][0] = 0; edgeNormals[2][1] = 0;
+    edgeNormals[3][0] = 0; edgeNormals[3][1] = 0;
+
+    forAll(edgeDirections, eI)
+    {
+        edgeDirections[eI] =
+            surface().points()()[edges[eI].end()]
+          - surface().points()()[edges[eI].start()];
+
+        normalDirections[eI] = labelList(2, 0);
+        for (label j = 0; j < 2; ++j)
+        {
+            const vector cross =
+                (faceNormals[edgeNormals[eI][j]] ^ edgeDirections[eI]);
+            const vector fC0tofE0 =
+                0.5*(max(surface().points()() + min(surface().points()())))
+              - surface().points()()[edges[eI].start()];
+
+            normalDirections[eI][j] =
+                (
+                    (
+                        (cross/(mag(cross) + VSMALL))
+                      & (fC0tofE0/(mag(fC0tofE0)+ VSMALL))
+                    )
+                  > 0.0
+                    ? 1
+                    : -1
+                );
+        }
+    }
+
+    labelListList featurePointNormals(4);
+    labelListList featurePointEdges(4);
+    forAll(featurePointNormals, pI)
+    {
+        labelList& ftPtEdges = featurePointEdges[pI];
+        ftPtEdges.setSize(2, 0);
+
+        label edgeI = 0;
+        forAll(edges, eI)
+        {
+            const edge& e = edges[eI];
+
+            if (e.start() == pI)
+            {
+                ftPtEdges[edgeI++] = eI;
+            }
+            else if (e.end() == pI)
+            {
+                ftPtEdges[edgeI++] = eI;
+            }
+        }
+
+        labelList& ftPtNormals = featurePointNormals[pI];
+        ftPtNormals.setSize(1, 0);
+
+        ftPtNormals[0] = edgeNormals[ftPtEdges[0]][0];
+    }
+
+    labelList regionEdges;
+
+    features.set
+    (
+        new extendedFeatureEdgeMesh
+        (
+            IOobject
+            (
+                surface().name(),
+                surface().instance(),
+                "extendedFeatureEdgeMesh",
+                surface().db(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            surface().points(),
+            edges,
+            4, 4, 4,
+            0, 0, 4, 4, // 4 flat edges
+            faceNormals,
+            List<extendedFeatureEdgeMesh::sideVolumeType>(4, mode_),
+            edgeDirections,
+            normalDirections,
+            edgeNormals,
+            featurePointNormals,
+            featurePointEdges,
+            regionEdges
+        )
+    );
+
+    return features;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.H
new file mode 100644
index 0000000000..3854a67d1f
--- /dev/null
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/searchablePlateFeatures.H
@@ -0,0 +1,120 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/>.
+
+Class
+    Foam::searchablePlateFeatures
+
+Description
+
+SourceFiles
+    searchablePlateFeatures.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchablePlateFeatures_H
+#define searchablePlateFeatures_H
+
+#include "searchableSurfaceFeatures.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class searchablePlateFeatures Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchablePlateFeatures
+:
+    public searchableSurfaceFeatures
+{
+private:
+
+        //- To initialise edges.
+        static edgeList calcEdges(const label[4][2]);
+
+
+    // Private Member Data
+
+        //- Which side of the plate to mesh
+        const extendedFeatureEdgeMesh::sideVolumeType mode_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        searchablePlateFeatures(const searchablePlateFeatures&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const searchablePlateFeatures&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("searchablePlateFeatures");
+
+
+    // Static data members
+
+        //- Edge to point addressing
+        static const edgeList edges;
+
+
+    // Constructors
+
+        //- Construct from searchable surface and dictionary
+        searchablePlateFeatures
+        (
+            const searchableSurface& surface,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~searchablePlateFeatures();
+
+
+    // Member Functions
+
+        //- Return true for a searchable plate having features
+        virtual bool hasFeatures() const
+        {
+            return true;
+        }
+
+        //- Return an extendedFeatureEdgeMesh containing the features
+        virtual autoPtr<extendedFeatureEdgeMesh> features() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.C
index 3af1dd80d8..c61b92f6f1 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.C
@@ -57,10 +57,19 @@ Foam::triSurfaceMeshFeatures::triSurfaceMeshFeatures
 )
 :
     searchableSurfaceFeatures(surface, dict),
-    includedAngle_(readScalar(dict.lookup("includedAngle")))
+    includedAngle_(readScalar(dict.lookup("includedAngle"))),
+    mode_
+    (
+        extendedFeatureEdgeMesh::sideVolumeTypeNames_
+        [
+            dict.lookupOrDefault<word>("meshableSide", "inside")
+        ]
+    )
 {
     Info<< indent
-        << "    Included angle = " << includedAngle_
+        << "    Included angle = " << includedAngle_ << nl
+        << "    Meshable region = "
+        << extendedFeatureEdgeMesh::sideVolumeTypeNames_[mode_]
         << endl;
 }
 
@@ -82,7 +91,12 @@ Foam::triSurfaceMeshFeatures::features() const
 
     surfaceFeatures sFeat(surfMesh, includedAngle_);
 
-    boolList surfBaffleRegions(surfMesh.patches().size(), false);
+    // @todo Need to read on a per region basis
+    boolList surfBaffleRegions
+    (
+        surfMesh.patches().size(),
+        (mode_ == extendedFeatureEdgeMesh::BOTH ? true : false)
+    );
 
     features.set
     (
@@ -90,8 +104,8 @@ Foam::triSurfaceMeshFeatures::features() const
         (
             sFeat,
             surface().db(),
-            surface().name() + ".extendedFeatureEdgeMesh"//,
-            //surfBaffleRegions
+            surface().name() + ".extendedFeatureEdgeMesh",
+            surfBaffleRegions
         )
     );
 
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.H b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.H
index 190b7859ce..3105a38943 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.H
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/searchableSurfaceFeatures/triSurfaceMeshFeatures.H
@@ -55,6 +55,9 @@ private:
 
         const scalar includedAngle_;
 
+        //- Which side of the surface to mesh
+        const extendedFeatureEdgeMesh::sideVolumeType mode_;
+
 
     // Private Member Functions
 
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
index debb7b0588..e5c928e298 100644
--- a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
+++ b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
@@ -78,11 +78,79 @@ Description
 #include "booleanSurface.H"
 #include "edgeIntersections.H"
 #include "meshTools.H"
+#include "labelPair.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+bool intersectSurfaces
+(
+    triSurface& surf,
+    edgeIntersections& edgeCuts
+)
+{
+    bool hasMoved = false;
+
+    for (label iter = 0; iter < 10; iter++)
+    {
+        Info<< "Determining intersections of surface edges with itself" << endl;
+
+        // Determine surface edge intersections. Allow surface to be moved.
+
+        // Number of iterations needed to resolve degenerates
+        label nIters = 0;
+        {
+            triSurfaceSearch querySurf(surf);
+
+            scalarField surfPointTol
+            (
+                1e-3*edgeIntersections::minEdgeLength(surf)
+            );
+
+            // Determine raw intersections
+            edgeCuts = edgeIntersections
+            (
+                surf,
+                querySurf,
+                surfPointTol
+            );
+
+            // Shuffle a bit to resolve degenerate edge-face hits
+            {
+                pointField points(surf.points());
+
+                nIters =
+                    edgeCuts.removeDegenerates
+                    (
+                        5,              // max iterations
+                        surf,
+                        querySurf,
+                        surfPointTol,
+                        points         // work array
+                    );
+
+                if (nIters != 0)
+                {
+                    // Update geometric quantities
+                    surf.movePoints(points);
+                    hasMoved = true;
+                }
+            }
+        }
+    }
+
+    if (hasMoved)
+    {
+        fileName newFile("surf.obj");
+        Info<< "Surface has been moved. Writing to " << newFile << endl;
+        surf.write(newFile);
+    }
+
+    return hasMoved;
+}
+
+
 // Keep on shuffling surface points until no more degenerate intersections.
 // Moves both surfaces and updates set of edge cuts.
 bool intersectSurfaces
@@ -238,6 +306,106 @@ label calcNormalDirection
 }
 
 
+void calcEdgeCuts
+(
+    triSurface& surf1,
+    triSurface& surf2,
+    const bool perturb,
+    edgeIntersections& edge1Cuts,
+    edgeIntersections& edge2Cuts
+)
+{
+    if (perturb)
+    {
+        intersectSurfaces
+        (
+            surf1,
+            edge1Cuts,
+            surf2,
+            edge2Cuts
+        );
+    }
+    else
+    {
+        triSurfaceSearch querySurf2(surf2);
+
+        Info<< "Determining intersections of surf1 edges with surf2 faces"
+            << endl;
+
+        edge1Cuts = edgeIntersections
+        (
+            surf1,
+            querySurf2,
+            1e-3*edgeIntersections::minEdgeLength(surf1)
+        );
+
+        triSurfaceSearch querySurf1(surf1);
+
+        Info<< "Determining intersections of surf2 edges with surf1 faces"
+            << endl;
+
+        edge2Cuts = edgeIntersections
+        (
+            surf2,
+            querySurf1,
+            1e-3*edgeIntersections::minEdgeLength(surf2)
+        );
+    }
+}
+
+
+void calcFeaturePoints(const pointField& points, const edgeList& edges)
+{
+    edgeMesh eMesh(points, edges);
+
+    const labelListList& pointEdges = eMesh.pointEdges();
+
+
+    // Get total number of feature points
+    label nFeaturePoints = 0;
+    forAll(pointEdges, pI)
+    {
+        const labelList& pEdges = pointEdges[pI];
+
+        if (pEdges.size() == 1)
+        {
+            nFeaturePoints++;
+        }
+    }
+
+
+    // Calculate addressing from feature point to cut point and cut edge
+    labelList featurePointToCutPoint(nFeaturePoints);
+    labelList featurePointToCutEdge(nFeaturePoints);
+
+    label nFeatPts = 0;
+    forAll(pointEdges, pI)
+    {
+        const labelList& pEdges = pointEdges[pI];
+
+        if (pEdges.size() == 1)
+        {
+            featurePointToCutPoint[nFeatPts] = pI;
+            featurePointToCutEdge[nFeatPts] = pEdges[0];
+            nFeatPts++;
+        }
+    }
+
+
+
+    label concaveStart = 0;
+    label mixedStart = 0;
+    label nonFeatureStart = nFeaturePoints;
+
+
+    labelListList featurePointNormals(nFeaturePoints);
+    labelListList featurePointEdges(nFeaturePoints);
+    labelList regionEdges;
+
+
+}
+
+
 int main(int argc, char *argv[])
 {
     argList::noParallel();
@@ -319,52 +487,33 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
-    if (args.optionFound("perturb"))
-    {
-        intersectSurfaces
-        (
-            surf1,
-            edge1Cuts,
-            surf2,
-            edge2Cuts
-        );
-    }
-    else
-    {
-        triSurfaceSearch querySurf2(surf2);
+    // Calculate the points where the edges are cut by the other surface
+    calcEdgeCuts
+    (
+        surf1,
+        surf2,
+        args.optionFound("perturb"),
+        edge1Cuts,
+        edge2Cuts
+    );
 
-        Info<< "Determining intersections of surf1 edges with surf2 faces"
-            << endl;
+    // Determine intersection edges from the edge cuts
+    surfaceIntersection inter
+    (
+        surf1,
+        edge1Cuts,
+        surf2,
+        edge2Cuts
+    );
 
-        edge1Cuts = edgeIntersections
-        (
-            surf1,
-            querySurf2,
-            1e-3*edgeIntersections::minEdgeLength(surf1)
-        );
-
-        triSurfaceSearch querySurf1(surf1);
-
-        Info<< "Determining intersections of surf2 edges with surf1 faces"
-            << endl;
-
-        edge2Cuts = edgeIntersections
-        (
-            surf2,
-            querySurf1,
-            1e-3*edgeIntersections::minEdgeLength(surf2)
-        );
-    }
-
-    // Determine intersection edges
-    surfaceIntersection inter(surf1, edge1Cuts, surf2, edge2Cuts);
-
-    fileName sFeatFileName =
+    fileName sFeatFileName
+    (
         surf1Name.lessExt().name()
       + "_"
       + surf2Name.lessExt().name()
       + "_"
-      + action;
+      + action
+    );
 
     label nFeatEds = inter.cutEdges().size();
 
@@ -393,6 +542,7 @@ int main(int argc, char *argv[])
         edgeDirections[cutEdgeI] = fE.vec(inter.cutPoints());
 
         normals.append(norm1);
+        eNormals.append(normals.size() - 1);
 
         if (surf1Baffle)
         {
@@ -416,9 +566,8 @@ int main(int argc, char *argv[])
             );
         }
 
-        eNormals.append(normals.size() - 1);
-
         normals.append(norm2);
+        eNormals.append(normals.size() - 1);
 
         if (surf2Baffle)
         {
@@ -443,8 +592,6 @@ int main(int argc, char *argv[])
             );
         }
 
-        eNormals.append(normals.size() - 1);
-
 
         if (surf1Baffle)
         {
@@ -509,6 +656,38 @@ int main(int argc, char *argv[])
 
 
     label internalStart = -1;
+    label nIntOrExt = 0;
+    label nFlat = 0;
+    label nOpen = 0;
+    label nMultiple = 0;
+
+    forAll(edgeNormals, eI)
+    {
+        label nEdNorms = edgeNormals[eI].size();
+
+        if (nEdNorms == 1)
+        {
+            nOpen++;
+        }
+        else if (nEdNorms == 2)
+        {
+            const vector& n0(normals[edgeNormals[eI][0]]);
+            const vector& n1(normals[edgeNormals[eI][1]]);
+
+            if ((n0 & n1) > extendedFeatureEdgeMesh::cosNormalAngleTol_)
+            {
+                nFlat++;
+            }
+            else
+            {
+                nIntOrExt++;
+            }
+        }
+        else if (nEdNorms > 2)
+        {
+            nMultiple++;
+        }
+    }
 
     if (validActions[action] == booleanSurface::UNION)
     {
@@ -520,7 +699,7 @@ int main(int argc, char *argv[])
         else
         {
             // All edges are external
-            internalStart = nFeatEds;
+            internalStart = nIntOrExt;
         }
     }
     else if (validActions[action] == booleanSurface::INTERSECTION)
@@ -528,7 +707,7 @@ int main(int argc, char *argv[])
         if (!invertedSpace)
         {
             // All edges are external
-            internalStart = nFeatEds;
+            internalStart = nIntOrExt;
         }
         else
         {
@@ -539,7 +718,7 @@ int main(int argc, char *argv[])
     else if (validActions[action] == booleanSurface::DIFFERENCE)
     {
         // All edges are external
-        internalStart = nFeatEds;
+        internalStart = nIntOrExt;
     }
     else
     {
@@ -569,6 +748,8 @@ int main(int argc, char *argv[])
         normalDirectionsTmp[i] = normalDirections[i];
     }
 
+    calcFeaturePoints(inter.cutPoints(), inter.cutEdges());
+
     extendedFeatureEdgeMesh feMesh
     (
         IOobject
@@ -582,18 +763,22 @@ int main(int argc, char *argv[])
         ),
         inter.cutPoints(),
         inter.cutEdges(),
+
         0,                  // concaveStart,
         0,                  // mixedStart,
         0,                  // nonFeatureStart,
+
         internalStart,      // internalStart,
-        nFeatEds,           // flatStart,
-        nFeatEds,           // openStart,
-        nFeatEds,           // multipleStart,
+        nIntOrExt,           // flatStart,
+        nIntOrExt + nFlat,   // openStart,
+        nIntOrExt + nFlat + nOpen,   // multipleStart,
+
         normalsTmp,
         normalVolumeTypesTmp,
         edgeDirections,
         normalDirectionsTmp,
         edgeNormalsTmp,
+
         labelListList(0),   // featurePointNormals,
         labelListList(0),   // featurePointEdges,
         labelList(0)        // regionEdges
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 58f321297b..7758293480 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -925,7 +925,7 @@ void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
         << "        open edges         : "
         << (fem.multipleStart()- fem.openStart()) << nl
         << "        multiply connected : "
-        << (fem.edges().size()- fem.multipleStart()) << nl;
+        << (fem.edges().size()- fem.multipleStart()) << endl;
 }
 
 
@@ -1223,17 +1223,34 @@ int main(int argc, char *argv[])
             << "        internal edges : " << newSet.nInternalEdges() << nl
             << endl;
 
-        //if (writeObj)
-        //{
-        //    newSet.writeObj("final");
-        //}
+        if (writeObj)
+        {
+            newSet.writeObj("final");
+        }
+
+        boolList surfBaffleRegions(surf.patches().size(), false);
+
+        wordList surfBaffleNames;
+        surfaceDict.readIfPresent("baffles", surfBaffleNames);
+
+        forAll(surf.patches(), pI)
+        {
+            const word& name = surf.patches()[pI].name();
+
+            if (findIndex(surfBaffleNames, name) != -1)
+            {
+                Info<< "Adding baffle region " << name << endl;
+                surfBaffleRegions[pI] = true;
+            }
+        }
 
         // Extracting and writing a extendedFeatureEdgeMesh
         extendedFeatureEdgeMesh feMesh
         (
             newSet,
             runTime,
-            sFeatFileName + ".extendedFeatureEdgeMesh"
+            sFeatFileName + ".extendedFeatureEdgeMesh",
+            surfBaffleRegions
         );
 
 
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
index 54ad050fc6..0d03912cee 100644
--- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
@@ -263,7 +263,8 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
 (
     const surfaceFeatures& sFeat,
     const objectRegistry& obr,
-    const fileName& sFeatFileName
+    const fileName& sFeatFileName,
+    const boolList& surfBaffleRegions
 )
 :
     regIOobject
@@ -337,7 +338,12 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
             // Check to see if the points have been already used
             if (faceMap[eFI] == -1)
             {
-                normalVolumeTypes_[nAdded++] = INSIDE;
+                normalVolumeTypes_[nAdded++] =
+                    (
+                        surfBaffleRegions[surf[eFI].region()]
+                      ? BOTH
+                      : INSIDE
+                    );
 
                 faceMap[eFI] = nAdded - 1;
             }
@@ -492,7 +498,7 @@ Foam::extendedFeatureEdgeMesh::classifyEdge
     const List<vector>& norms,
     const labelList& edNorms,
     const vector& fC0tofC1
-) const
+)
 {
     label nEdNorms = edNorms.size();
 
@@ -502,8 +508,8 @@ Foam::extendedFeatureEdgeMesh::classifyEdge
     }
     else if (nEdNorms == 2)
     {
-        const vector n0(norms[edNorms[0]]);
-        const vector n1(norms[edNorms[1]]);
+        const vector& n0(norms[edNorms[0]]);
+        const vector& n1(norms[edNorms[1]]);
 
         if ((n0 & n1) > cosNormalAngleTol_)
         {
@@ -1342,6 +1348,22 @@ void Foam::extendedFeatureEdgeMesh::writeObj
         meshTools::writeOBJ(regionStr, points()[e[1]]); verti++;
         regionStr << "l " << verti-1 << ' ' << verti << endl;
     }
+
+    OFstream edgeDirsStr(prefix + "_edgeDirections.obj");
+    Info<< "Writing edge directions to " << edgeDirsStr.name() << endl;
+
+    forAll(edgeDirections_, i)
+    {
+        const vector& eVec = edgeDirections_[i];
+        const edge& e = edges()[i];
+
+        meshTools::writeOBJ
+        (
+            edgeDirsStr,
+            points()[e.start()],
+            eVec + points()[e.start()]
+        );
+    }
 }
 
 
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
index dec5c41696..3695a39012 100644
--- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.H
@@ -119,14 +119,14 @@ public:
 
     static const Foam::NamedEnum<sideVolumeType, 4> sideVolumeTypeNames_;
 
+    //- Angular closeness tolerance for treating normals as the same
+    static scalar cosNormalAngleTol_;
+
 
 private:
 
     // Static data
 
-        //- Angular closeness tolerance for treating normals as the same
-        static scalar cosNormalAngleTol_;
-
         //- Index of the start of the convex feature points - static as 0
         static label convexStart_;
 
@@ -202,15 +202,6 @@ private:
         //  data for edges and normals.
         pointStatus classifyFeaturePoint(label ptI) const;
 
-        //- Classify the type of feature edge.  Requires face centre 0 to face
-        //  centre 1 vector to distinguish internal from external
-        edgeStatus classifyEdge
-        (
-            const List<vector>& norms,
-            const labelList& edNorms,
-            const vector& fC0tofC1
-        ) const;
-
         template<class Patch>
         void sortPointsAndEdges
         (
@@ -259,7 +250,8 @@ public:
         (
             const surfaceFeatures& sFeat,
             const objectRegistry& obr,
-            const fileName& sFeatFileName
+            const fileName& sFeatFileName,
+            const boolList& surfBaffleRegions
         );
 
         //- Construct from PrimitivePatch
@@ -434,6 +426,9 @@ public:
             //- Return the edgeStatus of a specified edge
             inline edgeStatus getEdgeStatus(label edgeI) const;
 
+            //- Return the baffle faces of a specified edge
+            inline PackedList<2> edgeBaffles(label edgeI) const;
+
             //- Demand driven construction of octree for feature points
             const indexedOctree<treeDataPoint>& pointTree() const;
 
@@ -468,6 +463,16 @@ public:
 
             friend Istream& operator>>(Istream& is, sideVolumeType& vt);
             friend Ostream& operator<<(Ostream& os, const sideVolumeType& vt);
+
+
+        //- Classify the type of feature edge.  Requires face centre 0 to face
+        //  centre 1 vector to distinguish internal from external
+        static edgeStatus classifyEdge
+        (
+            const List<vector>& norms,
+            const labelList& edNorms,
+            const vector& fC0tofC1
+        );
 };
 
 
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
index a0ba24daae..9259b06cb5 100644
--- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMeshI.H
@@ -264,4 +264,27 @@ Foam::extendedFeatureEdgeMesh::getEdgeStatus(label edgeI) const
 }
 
 
+inline Foam::PackedList<2> Foam::extendedFeatureEdgeMesh::edgeBaffles
+(
+    label edgeI
+) const
+{
+    const labelList& eNormals = edgeNormals_[edgeI];
+
+    DynamicList<label> edgeBaffles(eNormals.size());
+
+    forAll(eNormals, enI)
+    {
+        const label normI = eNormals[enI];
+
+        if (normalVolumeTypes_[normI])
+        {
+            edgeBaffles.append(normI);
+        }
+    }
+
+    return PackedList<2>(edgeBaffles);
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
index a233af491d..441aa3be52 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/edgeIntersections.C
@@ -117,19 +117,16 @@ void Foam::edgeIntersections::intersectEdges
             << " triangles ..." << endl;
     }
 
-    // Construct octree.
-    const indexedOctree<treeDataTriSurface>& tree = querySurf2.tree();
-
-
-    label nHits = 0;
-
+    pointField start(edgeLabels.size());
+    pointField end(edgeLabels.size());
+    vectorField edgeDirs(edgeLabels.size());
 
     // Go through all edges, calculate intersections
     forAll(edgeLabels, i)
     {
         label edgeI = edgeLabels[i];
 
-        if (debug && (i % 1000 == 0))
+        if (debug)// && (i % 1000 == 0))
         {
             Pout<< "Intersecting edge " << edgeI << " with surface" << endl;
         }
@@ -140,85 +137,81 @@ void Foam::edgeIntersections::intersectEdges
         const point& pEnd = points1[meshPoints[e.end()]];
 
         const vector eVec(pEnd - pStart);
-        const scalar eMag = mag(eVec);
-        const vector n(eVec/(eMag + VSMALL));
+        const vector n(eVec/(mag(eVec) + VSMALL));
 
-        // Smallish length for intersection calculation errors.
-        const point tolVec = 1e-6*eVec;
-
-        // Start tracking somewhat before pStart and upto somewhat after p1.
+        // Start tracking somewhat before pStart and up to somewhat after p1.
         // Note that tolerances here are smaller than those used to classify
         // hit below.
         // This will cause this hit to be marked as degenerate and resolved
         // later on.
-        point p0 = pStart - 0.5*surf1PointTol[e[0]]*n;
-        const point p1 = pEnd + 0.5*surf1PointTol[e[1]]*n;
-        const scalar maxS = mag(p1 - pStart);
+        start[i] = pStart - 0.5*surf1PointTol[e[0]]*n;
+        end[i] = pEnd + 0.5*surf1PointTol[e[1]]*n;
 
-        // Get all intersections of the edge with the surface
+        edgeDirs[i] = n;
+    }
 
-        DynamicList<pointIndexHit> currentIntersections(100);
-        DynamicList<label> currentIntersectionTypes(100);
+    List<List<pointIndexHit> > edgeIntersections;
+    querySurf2.findLineAll
+    (
+        start,
+        end,
+        edgeIntersections
+    );
 
-        while (true)
+    label nHits = 0;
+
+    // Classify the hits
+    forAll(edgeLabels, i)
+    {
+        const label edgeI = edgeLabels[i];
+
+        labelList& intersectionTypes = classification_[edgeI];
+        intersectionTypes.setSize(edgeIntersections[i].size(), -1);
+
+        this->operator[](edgeI).transfer(edgeIntersections[i]);
+
+        forAll(intersectionTypes, hitI)
         {
-            pointIndexHit pHit = tree.findLine(p0, p1);
+            const pointIndexHit& pHit = this->operator[](edgeI)[hitI];
+            label& hitType = intersectionTypes[hitI];
 
-            if (pHit.hit())
+            if (!pHit.hit())
             {
-                nHits++;
-
-                currentIntersections.append(pHit);
-
-                // Classify point on surface1 edge.
-                label edgeEnd = -1;
-
-                if (mag(pHit.hitPoint() - pStart) < surf1PointTol[e[0]])
-                {
-                    edgeEnd = 0;
-                }
-                else if (mag(pHit.hitPoint() - pEnd) < surf1PointTol[e[1]])
-                {
-                    edgeEnd = 1;
-                }
-                else if (mag(n & normals2[pHit.index()]) < alignedCos_)
-                {
-                    Pout<< "Flat angle edge:" << edgeI
-                        << " face:" << pHit.index()
-                        << " cos:" << mag(n & normals2[pHit.index()])
-                        << endl;
-                    edgeEnd = 2;
-                }
-
-                currentIntersectionTypes.append(edgeEnd);
-
-                if (edgeEnd == 1)
-                {
-                    // Close to end
-                    break;
-                }
-                else
-                {
-                    // Continue tracking. Shift by a small amount.
-                    p0 = pHit.hitPoint() + tolVec;
-
-                    if (((p0-pStart) & n) >= maxS)
-                    {
-                        break;
-                    }
-                }
+                continue;
             }
-            else
+
+            const edge& e = surf1.edges()[edgeI];
+
+            // Classify point on surface1 edge.
+            if (mag(pHit.hitPoint() - start[i]) < surf1PointTol[e[0]])
             {
-                // No hit.
-                break;
+                // Intersection is close to edge start
+                hitType = 0;
             }
+            else if (mag(pHit.hitPoint() - end[i]) < surf1PointTol[e[1]])
+            {
+                // Intersection is close to edge end
+                hitType = 1;
+            }
+            else if (mag(edgeDirs[i] & normals2[pHit.index()]) < alignedCos_)
+            {
+                // Edge is almost coplanar with the face
+
+                Pout<< "Flat angle edge:" << edgeI
+                    << " face:" << pHit.index()
+                    << " cos:" << mag(edgeDirs[i] & normals2[pHit.index()])
+                    << endl;
+
+                hitType = 2;
+            }
+
+            if (debug)
+            {
+                Info<< "    hit " << pHit << " classify = " << hitType << endl;
+            }
+
+            nHits++;
         }
-
-
-        // Done current edge. Transfer all data into *this
-        operator[](edgeI).transfer(currentIntersections);
-        classification_[edgeI].transfer(currentIntersectionTypes);
     }
 
     if (debug)
@@ -226,7 +219,6 @@ void Foam::edgeIntersections::intersectEdges
         Pout<< "Found " << nHits << " intersections of edges with surface ..."
             << endl;
     }
-
 }
 
 
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
index b54ad7143d..252bcdfbdc 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
@@ -202,7 +202,6 @@ void Foam::surfaceIntersection::storeIntersection
     DynamicList<point>& allCutPoints
 )
 {
-
     forAll(facesA, facesAI)
     {
         label faceA = facesA[facesAI];
@@ -1147,6 +1146,12 @@ const Foam::edgeList& Foam::surfaceIntersection::cutEdges() const
 }
 
 
+const Foam::labelPairLookup& Foam::surfaceIntersection::facePairToVertex() const
+{
+    return facePairToVertex_;
+}
+
+
 const Foam::labelPairLookup& Foam::surfaceIntersection::facePairToEdge() const
 {
     return facePairToEdge_;
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
index 68399c7e41..73638b4332 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
@@ -273,7 +273,7 @@ public:
 
         const edgeList& cutEdges() const;
 
-        //const labelPairLookup& facePairToVertex() const;
+        const labelPairLookup& facePairToVertex() const;
 
         const labelPairLookup& facePairToEdge() const;
 
diff --git a/tutorials/mesh/foamyHexMesh/mixerVessel/system/foamyHexMeshDict b/tutorials/mesh/foamyHexMesh/mixerVessel/system/foamyHexMeshDict
index 9ec0594223..dfbab846cb 100644
--- a/tutorials/mesh/foamyHexMesh/mixerVessel/system/foamyHexMeshDict
+++ b/tutorials/mesh/foamyHexMesh/mixerVessel/system/foamyHexMeshDict
@@ -59,9 +59,9 @@ surfaceConformation
             extendedFeatureEdgeMesh "vessel.extendedFeatureEdgeMesh";
             regions
             {
+                patch0  {}
                 patch1  {}
                 patch2  {}
-                patch3  {}
             }
         }
 
@@ -118,7 +118,7 @@ motionControl
 
     minimumCellSizeCoeff            0.1;
 
-    maxRefinementIterations         0;
+    maxRefinementIterations         1;
 
     maxSmoothingIterations          100;
 
@@ -142,7 +142,7 @@ motionControl
 
             regions
             {
-                patch1
+                patch0
                 {
                     surfaceCellSizeFunction     uniformValue;
                     uniformValueCoeffs
@@ -154,7 +154,7 @@ motionControl
                     uniformCoeffs{}
                 }
 
-                patch2
+                patch1
                 {
                     priority                    2;
                     surfaceCellSizeFunction     uniformValue;
@@ -170,7 +170,7 @@ motionControl
                     }
                 }
 
-                patch3
+                patch2
                 {
                     priority                    2;
                     surfaceCellSizeFunction     uniformValue;