diff --git a/applications/test/surfaceIntersection/Make/options b/applications/test/surfaceIntersection/Make/options
index 6ae6c04df1..b7fa8b2326 100644
--- a/applications/test/surfaceIntersection/Make/options
+++ b/applications/test/surfaceIntersection/Make/options
@@ -1,9 +1,7 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/edgeMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/triSurface/lnInclude
 
 EXE_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools -ledgeMesh
+    -lmeshTools
diff --git a/applications/test/surfaceIntersection/Test-surfaceIntersection.C b/applications/test/surfaceIntersection/Test-surfaceIntersection.C
index 8c0f83f641..ddf122f60a 100644
--- a/applications/test/surfaceIntersection/Test-surfaceIntersection.C
+++ b/applications/test/surfaceIntersection/Test-surfaceIntersection.C
@@ -34,7 +34,7 @@ Description
 #include "triSurface.H"
 #include "triSurfaceMesh.H"
 #include "surfaceIntersection.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 
 using namespace Foam;
 
@@ -172,7 +172,7 @@ int main(int argc, char *argv[])
     {
         Info<< "surf1-cuts: " << cuts.surf1EdgeCuts() << nl
             << "surf2-cuts: " << cuts.surf2EdgeCuts() << nl
-            << "face-pairs: " << cuts.facePairToEdge() << nl
+            << "face-pairs: " << cuts.facePairToEdgeId() << nl
             << "edges: " << cuts.cutEdges() << nl;
     }
 
@@ -198,7 +198,7 @@ int main(int argc, char *argv[])
         {
             Info<< "surf1-cuts: " << cuts.surf1EdgeCuts() << nl
                 << "surf2-cuts: " << cuts.surf2EdgeCuts() << nl
-                << "face-pairs: " << cuts.facePairToEdge() << nl
+                << "face-pairs: " << cuts.facePairToEdgeId() << nl
                 << "edges: " << cuts.cutEdges() << nl;
         }
     }
@@ -209,20 +209,7 @@ int main(int argc, char *argv[])
     if (points.size() || edges.size())
     {
         Info<<"write to " << outputFile << nl;
-
-        OFstream os(outputFile);
-
-        forAll(points, pointi)
-        {
-            const point& pt = points[pointi];
-            os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
-        }
-
-        forAll(edges, edgei)
-        {
-            const edge& e = edges[edgei];
-            os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
-        }
+        OBJstream(outputFile).write(edges, points);
     }
 
     Info<< "End\n" << endl;
diff --git a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
index 02d0b0c8e6..1d0dfae584 100644
--- a/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
+++ b/applications/utilities/surface/surfaceBooleanFeatures/surfaceBooleanFeatures.C
@@ -99,7 +99,7 @@ typedef CGAL::AABB_face_graph_triangle_primitive
 typedef CGAL::AABB_traits<K, Primitive> Traits;
 typedef CGAL::AABB_tree<Traits> Tree;
 
-typedef boost::optional<Tree::Intersection_and_primitive_id<Segment>::Type >
+typedef boost::optional<Tree::Intersection_and_primitive_id<Segment>::Type>
 Segment_intersection;
 
 #endif // NO_CGAL
@@ -477,7 +477,6 @@ label dupNonManifoldPoints(triSurface& s, labelList& pointMap)
     List<labelledTri> newFaces(s);
     label nNonManifold = 0;
 
-
     forAll(pf, pointI)
     {
         const labelList& pFaces = pf[pointI];
@@ -1257,10 +1256,10 @@ autoPtr<extendedFeatureEdgeMesh> createEdgeMesh
     const triSurface& s1 = surf1;
     const triSurface& s2 = surf2;
 
-    forAllConstIter(labelPairLookup, inter.facePairToEdge(), iter)
+    forAllConstIters(inter.facePairToEdgeId(), iter)
     {
-        const label& cutEdgeI = iter();
         const labelPair& facePair = iter.key();
+        const label cutEdgeI = iter.object();
 
         const edge& fE = inter.cutEdges()[cutEdgeI];
 
diff --git a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/surfaceFeaturesExtraction.C b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/surfaceFeaturesExtraction.C
index 95a1c1b028..28bd8f092d 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/surfaceFeaturesExtraction.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/surfaceFeaturesExtraction.C
@@ -73,7 +73,7 @@ Foam::surfaceFeaturesExtraction::method::New
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(methodName);
 
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    if (!cstrIter.found())
     {
         FatalIOErrorInFunction
         (
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 1ea773f16b..29187c7cfe 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -84,6 +84,8 @@ int main(int argc, char *argv[])
 
     forAllConstIter(dictionary, dict, iter)
     {
+        const word& dictName = iter().keyword();
+
         if (!iter().isDict())
         {
             continue;
@@ -106,7 +108,7 @@ int main(int argc, char *argv[])
         const word outputName =
             fileName
             (
-                surfaceDict.lookupOrDefault<word>("output", iter().keyword())
+                surfaceDict.lookupOrDefault<word>("output", dictName)
             ).lessExt();
 
         // The "surfaces" entry is normally optional, but if the sub-dictionary
@@ -115,7 +117,7 @@ int main(int argc, char *argv[])
         // additional switch.
         if
         (
-            iter().keyword() == "surfaces"  // mandatory
+            dictName == "surfaces"  // mandatory
          || surfaceDict.found("surfaces")   // or optional
         )
         {
@@ -123,14 +125,14 @@ int main(int argc, char *argv[])
         }
         else
         {
-            loader.select(iter().keyword());
+            loader.select(dictName);
         }
 
         if (loader.selected().empty())
         {
             FatalErrorInFunction
                 << "No surfaces specified/found for entry: "
-                << iter().keyword() << exit(FatalError);
+                << dictName << exit(FatalError);
         }
         // DebugVar(loader.available());
         // DebugVar(outputName);
@@ -153,7 +155,7 @@ int main(int argc, char *argv[])
         {
             FatalErrorInFunction
                 << "Problem loading surface(s) for entry: "
-                << iter().keyword() << exit(FatalError);
+                << dictName << exit(FatalError);
         }
 
         triSurface surf = surfPtr();
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
index d04bbc7958..362e7f99ee 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
@@ -134,7 +134,7 @@ surface2.nas
 // - If other dictionaries contain a 'surfaces' entry,
 //   it will be taken for the input.
 //
-surfaces
+dummyName
 {
     extractionMethod    extractFromSurface;
 
@@ -169,4 +169,39 @@ surfaces
         writeObj        yes;
 }
 
+
+// Handle single or multiple surfaces
+//
+// - If the dictionary is named 'surfaces', it must also contain a 'surfaces'
+//   entry (wordRe list).
+//
+// - If other dictionaries contain a 'surfaces' entry,
+//   it will be taken for the input.
+//
+surfaces
+{
+    extractionMethod    extractFromNone;
+
+    surfaces            (surface1.stl surface2.nas);
+
+    // Base output name (optional)
+    // output              surfaces;
+
+    // Generate additional features from self-intersect
+    selfIntersection    true;
+
+    // Tolerance for surface intersections
+    tolerance           1e-3;
+
+    extractFromNoneCoeffs
+    {
+        includedAngle   0;
+    }
+
+    // Write options
+
+    // Write features to obj format for postprocessing
+    writeObj        yes;
+}
+
 // ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C b/applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C
index 65dd086cd3..47a83a09f2 100644
--- a/applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C
+++ b/applications/utilities/surface/surfaceSplitNonManifolds/surfaceSplitNonManifolds.C
@@ -129,8 +129,7 @@ void dumpFaces
     const Map<label>& connectedFaces
 )
 {
-    Info<< "Dumping connectedFaces as Lightwave .obj file to " << fName
-        << "\nThis can be visualized with e.g. javaview (www.javaview.de)\n\n";
+    Info<< "Dumping connectedFaces as .obj file to " << fName << nl;
 
     OFstream os(fName);
 
diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C b/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C
index 87405317a2..8491d6dee1 100644
--- a/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C
+++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.C
@@ -27,56 +27,37 @@ License
 #include "triSurface.H"
 #include "surfaceIntersection.H"
 #include "meshTools.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
 defineTypeNameAndDebug(edgeSurface, 0);
-}
 
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-// Write whole pointField and edges to stream
-void Foam::edgeSurface::writeOBJ
-(
-    const pointField& points,
-    const edgeList& edges,
-    Ostream& os
-)
+// file-scope
+// Write points in obj format
+static void writeObjPoints(const UList<point>& pts, Ostream& os)
 {
-    forAll(points, pointi)
+    forAll(pts, i)
     {
-        const point& pt = points[pointi];
-
-        os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
-    }
-    forAll(edges, edgeI)
-    {
-        const edge& e = edges[edgeI];
-
-        os << "l " << e.start()+1 << ' ' << e.end()+1 << endl;
+        const point& pt = pts[i];
+        os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
     }
 }
 
 
 // Write whole pointField and selected edges to stream
-void Foam::edgeSurface::writeOBJ
+void writeObjEdges
 (
-    const pointField& points,
+    const UList<point>& points,
     const edgeList& edges,
     const labelList& edgeLabels,
     Ostream& os
 )
 {
-    forAll(points, pointi)
-    {
-        const point& pt = points[pointi];
+    writeObjPoints(points, os);
 
-        os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
-    }
     forAll(edgeLabels, i)
     {
         const edge& e = edges[edgeLabels[i]];
@@ -85,6 +66,10 @@ void Foam::edgeSurface::writeOBJ
     }
 }
 
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 // Pointedges in edgeSurface indices only.
 void Foam::edgeSurface::calcPointEdges()
@@ -253,29 +238,16 @@ Foam::edgeSurface::edgeSurface
     }
 
 
-
-
     // Add intersection edges to faceEdges
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    forAllConstIter(labelPairLookup, inter.facePairToEdge(), iter)
+    forAllConstIters(inter.facePairToEdgeId(), iter)
     {
+        // The faceId from the correct surface
+        const label facei = iter.key()[isFirstSurface ? 0 : 1];
+
         // Edge label in intersection
-        const label edgeI = iter();
-
-        // Get the face from the correct surface
-        const FixedList<label, 2>& twoFaces = iter.key();
-
-        label facei;
-
-        if (isFirstSurface)
-        {
-            facei = twoFaces[0];
-        }
-        else
-        {
-            facei = twoFaces[1];
-        }
+        const label edgeI = iter.object();
 
         // Store on face-edge addressing. (note: offset edge)
         allFaceEdges[facei].append(edgeI + nSurfaceEdges_);
@@ -312,18 +284,17 @@ Foam::edgeSurface::edgeSurface
                     << " to " << faceFName << endl;
 
                 OFstream fStream(faceFName);
-                writeOBJ(points_, edges_, fEdges, fStream);
+                writeObjEdges(points_, edges_, fEdges, fStream);
             }
         }
 
         Pout<< "edgeSurface : Dumping edges to edges.obj" << endl;
-        OFstream eStream("edges.obj");
-        writeOBJ(points_, edges_, eStream);
+        OBJstream("edges.obj").write(edges_, points_);
 
         Pout<< "edgeSurface : Dumping intersectionEdges to"
             << " intersectionEdges.obj" << endl;
-        OFstream intEdgesStream("intersectionEdges.obj");
 
+        OFstream intEdgesStream("intersectionEdges.obj");
         labelList edgeLabels(edges_.size() - nSurfaceEdges_);
 
         label i = 0;
@@ -332,7 +303,7 @@ Foam::edgeSurface::edgeSurface
             edgeLabels[i++] = edgeI;
         }
 
-        writeOBJ(points_, edges_, edgeLabels, intEdgesStream);
+        writeObjEdges(points_, edges_, edgeLabels, intEdgesStream);
     }
 }
 
diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.H b/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.H
index d8564be949..76c83ee2ed 100644
--- a/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.H
+++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/edgeSurface.H
@@ -93,25 +93,12 @@ private:
         //- From face to our edges_
         labelListList faceEdges_;
 
-
         //- Constructed from above: pointEdges
         labelListList pointEdges_;
 
 
     // Private Member Functions
 
-        //- Dump edges in obj format
-        static void writeOBJ(const pointField&, const edgeList&, Ostream&);
-
-        //- Dump selected edges in obj format
-        static void writeOBJ
-        (
-            const pointField&,
-            const edgeList&,
-            const labelList&,
-            Ostream&
-        );
-
         //- Calculate pointEdges
         void calcPointEdges();
 
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
index 737387a5e1..5592d95362 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
@@ -25,7 +25,7 @@ License
 
 #include "surfaceIntersection.H"
 #include "triSurfaceSearch.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 #include "labelPairHashes.H"
 #include "triSurface.H"
 #include "pointIndexHit.H"
@@ -47,7 +47,7 @@ void Foam::surfaceIntersection::setOptions(const dictionary& dict)
 {
     dict.readIfPresent("tolerance",       tolerance_);
     dict.readIfPresent("allowEdgeHits",   allowEdgeHits_);
-    dict.readIfPresent("avoidDuplicates", avoidDuplicates_);
+    dict.readIfPresent("snap",            snapToEnd_);
     dict.readIfPresent("warnDegenerate",  warnDegenerate_);
 }
 
@@ -101,42 +101,64 @@ void Foam::surfaceIntersection::storeIntersection
             }
         }
 
-        labelPairLookup::const_iterator iter = facePairToVertex_.find(twoFaces);
 
-        if (iter == facePairToVertex_.end())
+        // Get existing edge, or create a null edge (with -1)
+        edge& thisEdge = facePairToEdge_(twoFaces);
+        const label pointCount = thisEdge.count();
+
+        if (pointCount == 0)
         {
-            // New intersection. Store face-face intersection.
+            // First intersection of the faces - record it.
+            thisEdge.insert(cutPointId);
+
             if (debug & 4)
             {
                 Pout<< "intersect faces " << twoFaces
                     << " point-1: " << cutPointId << " = "
                     << allCutPoints[cutPointId] << endl;
             }
-
-            facePairToVertex_.insert(twoFaces, cutPointId);
+            continue;
         }
-        else if (*iter == cutPointId)
+        else if (pointCount == 2)
         {
-            // Avoid creating an edge if cutPointId had already been used
-
+            // This occurs for ugly surfaces with shards that result in multiple
+            // cuts very near a snapped end point.
             if (debug & 4)
             {
-                Pout<< "intersect faces " << twoFaces
-                    << " dup-point: " << cutPointId << endl;
+                Pout<< "suppressed double intersection " << twoFaces
+                    << thisEdge << endl;
             }
+            continue;
         }
-        else
+
+        if (thisEdge.insert(cutPointId))
         {
-            const label nextEdgeId = allCutEdges.size();
-            const edge nextEdge(*iter, cutPointId, true);
+            // Second intersection of the faces - this is an edge,
+            // with special treatment:
+            // - avoid duplicate points: addressed by the insert() above
+            // - avoid degenerate lengths
+            // - avoid duplicate edges - can occur with really dirty geometry
 
-            // Second occurrence of surf1-surf2 intersection.
-            // Or rather the face on surf1 intersects a face on
-            // surface2 twice -> we found edge.
-
-            // Check whether perhaps degenerate
-            if (nextEdge.mag(allCutPoints) < SMALL)
+            if (edgeToId_.found(thisEdge))
             {
+                // Already have this edgeId, but not for this intersection.
+                thisEdge.sort();
+                if (facePairToEdge_.insert(twoFaces, thisEdge))
+                {
+                    if (debug & 4)
+                    {
+                        Pout<< "reuse edge - faces " << twoFaces << " edge#"
+                            << edgeToId_[thisEdge] << " edge " << thisEdge
+                            << " = " << thisEdge.line(allCutPoints)
+                            << endl;
+                    }
+                }
+            }
+            else if (thisEdge.mag(allCutPoints) < SMALL)
+            {
+                // Degenerate length
+                // - eg, end snapping was disabled or somehow failed.
+
                 // Don't normally emit warnings, since these also arise for
                 // manifold connections. For example,
                 //
@@ -147,10 +169,6 @@ void Foam::surfaceIntersection::storeIntersection
                 //
                 // The plane is correctly pierced at the '.' by both edge-1
                 // and edge-2, which belong to the same originating face.
-                //
-                // Unfortunately cannot suppress the second hit either, since
-                // it might already have been used for another face-pair
-                // intersection.
 
                 // Filter/merge away the extraneous points later.
                 if (warnDegenerate_ > 0)
@@ -159,31 +177,55 @@ void Foam::surfaceIntersection::storeIntersection
                     WarningInFunction
                         << "Degenerate edge between faces " << twoFaces
                         << " on 1st/2nd surface with points "
-                        << nextEdge.line(allCutPoints)
+                        << thisEdge.line(allCutPoints)
                         << endl;
                 }
                 else if (debug & 4)
                 {
                     Pout<< "degenerate edge face-pair " << twoFaces << " "
-                        << *iter << " point " << allCutPoints[*iter]
+                        << thisEdge[0] << " point " << allCutPoints[thisEdge[0]]
                         << endl;
                 }
+
+                // This is a failed edge - undo this second interaction
+                thisEdge.erase(cutPointId);
             }
-            else if (facePairToEdge_.insert(twoFaces, nextEdgeId))
+            else
             {
-                // Record complete (line) intersection of two faces
+                // This is a new edge.
+                const label edgeId = allCutEdges.size();
 
-                allCutEdges.append(nextEdge);
-
-                if (debug & 4)
+                if (facePairToEdgeId_.insert(twoFaces, edgeId))
                 {
-                    Pout<< "create edge - faces " << twoFaces << " edge#"
-                        << nextEdgeId << " edge " << nextEdge
-                        << " = " << nextEdge.line(allCutPoints)
-                        << endl;
+                    // Record complete (line) intersection of two faces
+                    thisEdge.sort();
+                    edgeToId_.insert(thisEdge, edgeId);
+                    allCutEdges.append(thisEdge);
+
+                    if (debug & 4)
+                    {
+                        Pout<< "create edge - faces " << twoFaces << " edge#"
+                            << edgeId << " edge " << thisEdge
+                            << " = " << thisEdge.line(allCutPoints)
+                            << endl;
+                    }
+                }
+                else
+                {
+                    // Faces already had an intersection
+                    // This should not fail, but for safety.
+                    Info<<"WARN " << twoFaces
+                        << " already intersected= " << thisEdge << endl;
+                    thisEdge.erase(cutPointId);
                 }
             }
         }
+        else
+        {
+            // Duplicate point - usually zero-length edge from snapping
+            // - can discard this face/face interaction entirely
+            facePairToEdge_.erase(twoFaces);
+        }
     }
 }
 
@@ -212,7 +254,7 @@ void Foam::surfaceIntersection::classifyHit
     List<DynamicList<label>>& surfEdgeCuts
 )
 {
-    const edge& e = surf1.edges()[edgeI];
+    const edge& e1 = surf1.edges()[edgeI];
 
     const labelList& facesA = surf1.edgeFaces()[edgeI];
 
@@ -233,10 +275,10 @@ void Foam::surfaceIntersection::classifyHit
     const label edgeEnd =
         classify
         (
-            surf1PointTol[e.start()],
-            surf1PointTol[e.end()],
+            surf1PointTol[e1.start()],
+            surf1PointTol[e1.end()],
             pHit.hitPoint(),
-            e,
+            e1,
             surf1Pts
         );
 
@@ -248,8 +290,8 @@ void Foam::surfaceIntersection::classifyHit
             if (debug & 2)
             {
                 Pout<< "hit-type[1] " << pHit.hitPoint() << " is surf1:"
-                    << " end point of edge[" << edgeI << "] " << e
-                    << "==" << e.line(surf1Pts)
+                    << " end point of edge[" << edgeI << "] " << e1
+                    << "==" << e1.line(surf1Pts)
                     << " surf2: vertex " << f2[nearLabel]
                     << " coord:" << surf2Pts[f2[nearLabel]]
                     << " - suppressed" << endl;
@@ -258,7 +300,6 @@ void Foam::surfaceIntersection::classifyHit
         else
         {
             // 2. Edge hits point. Cut edge with new point.
-            bool cached = false;
             label cutPointId = -1;
             const label nearVert = f2[nearLabel];
 
@@ -269,22 +310,25 @@ void Foam::surfaceIntersection::classifyHit
             {
                 const point& nearPt = surf1Pts[nearVert];
 
-                if (mag(pHit.hitPoint() - nearPt) < surf1PointTol[nearVert])
+                if
+                (
+                    mag(pHit.hitPoint() - nearPt)
+                  < surf1PointTol[nearVert]
+                )
                 {
                     cutPointId = allCutPoints.size();
 
-                    if (avoidDuplicates_)
+                    if (snapToEnd_)
                     {
-                        if (edgeEndAsCut_.insert(nearVert, cutPointId))
+                        if (snappedEnds_.insert(nearVert, cutPointId))
                         {
-                            // First time with this end-point
+                            // Initial snap
                             allCutPoints.append(nearPt);
                         }
                         else
                         {
-                            // Already seen this end point
-                            cutPointId = edgeEndAsCut_[nearVert];
-                            cached = true;
+                            // Already snapped this point.
+                            cutPointId = snappedEnds_[nearVert];
                         }
                     }
                     else
@@ -297,11 +341,11 @@ void Foam::surfaceIntersection::classifyHit
             if (debug & 2)
             {
                 Pout<< "hit-type[2] " << pHit.hitPoint() << " is surf1:"
-                    << " from edge[" << edgeI << "] " << e
+                    << " from edge[" << edgeI << "] " << e1
                     << " surf2: vertex " << f2[nearLabel]
                     << " coord:" << surf2Pts[f2[nearLabel]]
                     << " - "
-                    << (cached ? "cached" : "stored") << endl;
+                    << (cutPointId >= 0 ? "snapped" : "stored") << endl;
             }
 
             if (cutPointId == -1)
@@ -339,7 +383,7 @@ void Foam::surfaceIntersection::classifyHit
 
             const label edge2I = getEdge(surf2, surf2Facei, nearLabel);
             const edge& e2 = surf2.edges()[edge2I];
-            const label nearVert  = (edgeEnd == 0 ? e.start() : e.end());
+            const label nearVert = e1[edgeEnd];
 
             label cutPointId = -1;
 
@@ -367,22 +411,23 @@ void Foam::surfaceIntersection::classifyHit
 
                     if
                     (
-                        mag(pHit.hitPoint() - nearPt) < surf1PointTol[nearVert]
+                        mag(pHit.hitPoint() - nearPt)
+                      < surf1PointTol[nearVert]
                     )
                     {
                         cutPointId = allCutPoints.size();
 
-                        if (avoidDuplicates_)
+                        if (snapToEnd_)
                         {
-                            if (edgeEndAsCut_.insert(nearVert, cutPointId))
+                            if (snappedEnds_.insert(nearVert, cutPointId))
                             {
-                                // First time with this end-point
+                                // Initial snap
                                 allCutPoints.append(nearPt);
                             }
                             else
                             {
-                                // Already seen this end point
-                                cutPointId = edgeEndAsCut_[nearVert];
+                                // Already snapped this point.
+                                cutPointId = snappedEnds_[nearVert];
                                 handling = 2;  // cached
                             }
                         }
@@ -401,8 +446,8 @@ void Foam::surfaceIntersection::classifyHit
             if (debug & 2)
             {
                 Pout<< "hit-type[3] " << pHit.hitPoint() << " is surf1:"
-                    << " end point of edge[" << edgeI << "] " << e
-                    << "==" << e.line(surf1Pts)
+                    << " end point of edge[" << edgeI << "] " << e1
+                    << "==" << e1.line(surf1Pts)
                     << " surf2: edge[" << edge2I << "] " << e2
                     << " coords:" << e2.line(surf2Pts)
                     << " - "
@@ -482,9 +527,9 @@ void Foam::surfaceIntersection::classifyHit
                     if (edgeEdgeIntersection_.insert(intersect))
                     {
                         handling = 1;
-                        forAll(e, edgepti)
+                        forAll(e1, edgepti)
                         {
-                            const label endId = e[edgepti];
+                            const label endId = e1[edgepti];
                             const point& nearPt = surf1Pts[endId];
 
                             if
@@ -495,9 +540,9 @@ void Foam::surfaceIntersection::classifyHit
                             {
                                 cutPointId = allCutPoints.size();
 
-                                if (avoidDuplicates_)
+                                if (snapToEnd_)
                                 {
-                                    if (edgeEndAsCut_.insert(endId, cutPointId))
+                                    if (snappedEnds_.insert(endId, cutPointId))
                                     {
                                         // First time with this end-point
                                         allCutPoints.append(nearPt);
@@ -505,7 +550,7 @@ void Foam::surfaceIntersection::classifyHit
                                     else
                                     {
                                         // Already seen this end point
-                                        cutPointId = edgeEndAsCut_[endId];
+                                        cutPointId = snappedEnds_[endId];
                                         handling = 2;  // cached
                                     }
                                 }
@@ -524,8 +569,8 @@ void Foam::surfaceIntersection::classifyHit
             if (debug & 2)
             {
                 Pout<< "hit-type[4] " << pHit.hitPoint() << " is surf1:"
-                    << " from edge[" << edgeI << "] " << e
-                    << "==" << e.line(surf1Pts)
+                    << " from edge[" << edgeI << "] " << e1
+                    << "==" << e1.line(surf1Pts)
                     << " surf2: edge[" << edge2I << "] " << e2
                     << " coords:" << e2.line(surf2Pts)
                     << " - "
@@ -549,18 +594,28 @@ void Foam::surfaceIntersection::classifyHit
 
             if (handling)
             {
+                const vector eVec = e1.unitVec(surf1Pts);
+
                 const labelList& facesB = surf2.edgeFaces()[edge2I];
                 forAll(facesB, faceBI)
                 {
-                    storeIntersection
+                    // Intersecting edge should be non-coplanar with face
+                    if
                     (
-                        cutFrom,
-                        facesA,
-                        facesB[faceBI],
-                        allCutPoints,
-                        cutPointId,
-                        allCutEdges
-                    );
+                        mag((surf2.faceNormals()[facesB[faceBI]] & eVec))
+                      > 0.01
+                    )
+                    {
+                        storeIntersection
+                        (
+                            cutFrom,
+                            facesA,
+                            facesB[faceBI],
+                            allCutPoints,
+                            cutPointId,
+                            allCutEdges
+                        );
+                    }
                 }
             }
         }
@@ -577,8 +632,8 @@ void Foam::surfaceIntersection::classifyHit
 
             // Vertex on/near surf2; vertex away from surf2
             // otherVert on outside of surf2
-            const label nearVert  = (edgeEnd == 0 ? e.start() : e.end());
-            const label otherVert = (edgeEnd == 0 ? e.end() : e.start());
+            const label nearVert  = (edgeEnd == 0 ? e1.start() : e1.end());
+            const label otherVert = (edgeEnd == 0 ? e1.end() : e1.start());
 
             const point& nearPt  = surf1Pts[nearVert];
             const point& otherPt = surf1Pts[otherVert];
@@ -592,9 +647,9 @@ void Foam::surfaceIntersection::classifyHit
                 bool cached = false;
 
                 label cutPointId = allCutPoints.size();
-                if (avoidDuplicates_)
+                if (snapToEnd_)
                 {
-                    if (edgeEndAsCut_.insert(nearVert, cutPointId))
+                    if (snappedEnds_.insert(nearVert, cutPointId))
                     {
                         // First time with this end-point
                         allCutPoints.append(nearPt);
@@ -602,7 +657,7 @@ void Foam::surfaceIntersection::classifyHit
                     else
                     {
                         // Already seen this end point
-                        cutPointId = edgeEndAsCut_[nearVert];
+                        cutPointId = snappedEnds_[nearVert];
                         cached = true;
                     }
                 }
@@ -617,8 +672,8 @@ void Foam::surfaceIntersection::classifyHit
                 {
                     Pout<< "hit-type[5] " << pHit.hitPoint()
                         << " shifted to " << nearPt
-                        << " from edge[" << edgeI << "] " << e
-                        << "==" << e.line(surf1Pts)
+                        << " from edge[" << edgeI << "] " << e1
+                        << "==" << e1.line(surf1Pts)
                         << " hits surf2 face[" << surf2Facei << "]"
                         << " - "
                         << (cached ? "cached" : "stored") << endl;
@@ -640,7 +695,7 @@ void Foam::surfaceIntersection::classifyHit
                 if (debug & 2)
                 {
                     Pout<< "hit-type[5] " << pHit.hitPoint()
-                        << " from edge[" << edgeI << "] " << e
+                        << " from edge[" << edgeI << "] " << e1
                         << " hits inside of surf2 face[" << surf2Facei << "]"
                         << " - discarded" << endl;
                 }
@@ -652,8 +707,8 @@ void Foam::surfaceIntersection::classifyHit
             if (debug & 2)
             {
                 Pout<< "hit-type[6] " << pHit.hitPoint()
-                    << " from edge[" << edgeI << "] " << e
-                    << "==" << e.line(surf1Pts)
+                    << " from edge[" << edgeI << "] " << e1
+                    << "==" << e1.line(surf1Pts)
                     << " hits surf2 face[" << surf2Facei << "]"
                     << " - stored" << endl;
             }
@@ -680,8 +735,8 @@ void Foam::surfaceIntersection::classifyHit
 // Cut all edges of surf1 with surf2. Sets
 // - cutPoints          : coordinates of cutPoints
 // - cutEdges           : newly created edges between cutPoints
-// - facePairToVertex   : hash from face1I and face2I to (first) cutPoint
-// - facePairToEdge     : hash from face1I and face2I to cutEdge
+// - facePairToVertex   : hash from face1I and face2I to edge
+// - facePairToEdgeId   : hash from face1I and face2I to index in cutEdge
 // - surfEdgeCuts       : gives for each edge the cutPoints
 //                        (in order from start to end)
 //
@@ -716,7 +771,7 @@ void Foam::surfaceIntersection::doCutEdges
         // An edge may intersect multiple faces
         // - mask out faces that have already been hit before trying again
         // - never intersect with faces attached to the edge itself
-        DynamicList<label> maskFaces(8);
+        DynamicList<label> maskFaces(32);
 
         treeDataTriSurface::findAllIntersectOp
             allIntersectOp(searchTree, maskFaces);
@@ -732,8 +787,12 @@ void Foam::surfaceIntersection::doCutEdges
             const point ptEnd =
                 surf1Pts[e.end()]   + 0.5*surf1PointTol[e.end()]*edgeVec;
 
-            // Never intersect with faces attached to the edge itself
-            maskFaces = surf1.edgeFaces()[edgeI];
+            // Never intersect with faces attached directly to the edge itself,
+            // nor with faces attached to its end points. This mask contains
+            // some duplicates, but filtering them out is less efficient.
+            maskFaces = surf1.pointFaces()[e.start()];
+            maskFaces.append(surf1.pointFaces()[e.end()]);
+
             while (true)
             {
                 pointIndexHit pHit = searchTree.findLine
@@ -833,24 +892,110 @@ void Foam::surfaceIntersection::doCutEdges
 
     // These temporaries are now unneeded:
     edgeEdgeIntersection_.clear();
-    edgeEndAsCut_.clear();
+    snappedEnds_.clear();
 
     intersection::setPlanarTol(oldTol);
 }
 
 
+void Foam::surfaceIntersection::joinDisconnected
+(
+    DynamicList<edge>& allCutEdges
+)
+{
+    // This simple heuristic seems to work just as well (or better) than
+    // more complicated schemes
+    //
+    // For any face/face intersection that only appears once,
+    // consider which other faces/points are involved and connect between
+    // those points.
+    // Just do a simple connect-the-dots?
+
+    Pair<Map<labelPairHashSet>> missedFacePoint;
+
+    // Stage 1:
+    // - Extract "faceId -> (faceId, pointId)"
+    //   for all face/face pairs that only have one interaction
+    forAllConstIters(facePairToEdge_, iter)
+    {
+        const labelPair& twoFaces = iter.key();
+        const edge& e = iter.object();
+
+        if (e.count() == 1)
+        {
+            // minVertex = -1 (unused), maxVertex = pointId
+            const label pointId = e.maxVertex();
+
+            missedFacePoint[0](twoFaces[0]).insert
+            (
+                labelPair(twoFaces[1], pointId)
+            );
+
+            missedFacePoint[1](twoFaces[1]).insert
+            (
+                labelPair(twoFaces[0], pointId)
+            );
+        }
+    }
+
+
+    // Stage 2:
+    // - anything with two cross-interactions could cause a new edge:
+
+    edgeHashSet newEdges;
+    forAll(missedFacePoint, sidei)
+    {
+        const auto& mapping = missedFacePoint[sidei];
+
+        forAllConstIters(mapping, iter)
+        {
+            const auto& connect = iter.object();
+
+            if (connect.size() == 2)
+            {
+                // exactly two face/face cross-interactions
+
+                edge e;
+                for (const auto& facePoint : connect)
+                {
+                    e.insert(facePoint.second());
+                }
+                e.sort();
+
+                // Only consider edges with two unique ends,
+                // and do not introduce duplicates
+                if (e.count() == 2 && !edgeToId_.found(e))
+                {
+                    newEdges.insert(e);
+                }
+            }
+        }
+    }
+
+    label edgeId = allCutEdges.size();
+    edgeList newEdgesLst = newEdges.sortedToc();
+    for (const auto& e : newEdgesLst)
+    {
+        // Record complete (line) intersection of two faces
+        allCutEdges.append(e);
+        edgeToId_.insert(e, edgeId);
+        ++edgeId;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::surfaceIntersection::surfaceIntersection()
 :
     tolerance_(1e-3),
     allowEdgeHits_(true),
-    avoidDuplicates_(true),
+    snapToEnd_(true),
     warnDegenerate_(0),
     cutPoints_(0),
     cutEdges_(0),
-    facePairToVertex_(0),
     facePairToEdge_(0),
+    facePairToEdgeId_(0),
     surf1EdgeCuts_(0),
     surf2EdgeCuts_(0)
 {}
@@ -865,12 +1010,12 @@ Foam::surfaceIntersection::surfaceIntersection
 :
     tolerance_(1e-3),
     allowEdgeHits_(true),
-    avoidDuplicates_(true),
+    snapToEnd_(true),
     warnDegenerate_(0),
     cutPoints_(0),
     cutEdges_(0),
-    facePairToVertex_(2*max(query1.surface().size(), query2.surface().size())),
     facePairToEdge_(2*max(query1.surface().size(), query2.surface().size())),
+    facePairToEdgeId_(2*max(query1.surface().size(), query2.surface().size())),
     surf1EdgeCuts_(0),
     surf2EdgeCuts_(0)
 {
@@ -931,6 +1076,9 @@ Foam::surfaceIntersection::surfaceIntersection
         edgeCuts2
     );
 
+    // join disconnected intersection points
+    joinDisconnected(allCutEdges);
+
     // Transfer to straight label(List)List
     transfer(edgeCuts2, surf2EdgeCuts_);
     cutEdges_.transfer(allCutEdges);
@@ -946,8 +1094,7 @@ Foam::surfaceIntersection::surfaceIntersection
         Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
             << endl;
 
-        OFstream intStream("intEdges.obj");
-        writeOBJ(cutPoints_, cutEdges_, intStream);
+        OBJstream("intEdges.obj").write(cutEdges_, cutPoints_);
 
         // Dump all cut edges to files
         Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
@@ -960,9 +1107,9 @@ Foam::surfaceIntersection::surfaceIntersection
     }
 
     // Temporaries
-    facePairToVertex_.clear();
+    facePairToEdge_.clear();
 
-    // // Cleanup any duplicate cuts?
+    // Cleanup any duplicate cuts?
     // mergeEdges();
 }
 
@@ -975,12 +1122,12 @@ Foam::surfaceIntersection::surfaceIntersection
 :
     tolerance_(1e-3),
     allowEdgeHits_(true),
-    avoidDuplicates_(true),
+    snapToEnd_(true),
     warnDegenerate_(0),
     cutPoints_(0),
     cutEdges_(0),
-    facePairToVertex_(2*query1.surface().size()),
     facePairToEdge_(2*query1.surface().size()),
+    facePairToEdgeId_(2*query1.surface().size()),
     surf1EdgeCuts_(0),
     surf2EdgeCuts_(0)
 {
@@ -1013,6 +1160,9 @@ Foam::surfaceIntersection::surfaceIntersection
         edgeCuts1
     );
 
+    // join disconnected intersection points
+    joinDisconnected(allCutEdges);
+
     // Transfer to straight label(List)List
     transfer(edgeCuts1, surf1EdgeCuts_);
     cutEdges_.transfer(allCutEdges);
@@ -1039,8 +1189,7 @@ Foam::surfaceIntersection::surfaceIntersection
         Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
             << endl;
 
-        OFstream intStream("intEdges.obj");
-        writeOBJ(cutPoints_, cutEdges_, intStream);
+        OBJstream("intEdges.obj").write(cutEdges_, cutPoints_);
 
         // Dump all cut edges to files
         Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
@@ -1049,7 +1198,7 @@ Foam::surfaceIntersection::surfaceIntersection
     }
 
     // Temporaries
-    facePairToVertex_.clear();
+    facePairToEdge_.clear();
 
     // // Cleanup any duplicate cuts?
     // mergeEdges();
@@ -1066,12 +1215,12 @@ Foam::surfaceIntersection::surfaceIntersection
 :
     tolerance_(1e-3),
     allowEdgeHits_(true),
-    avoidDuplicates_(true),
+    snapToEnd_(true),
     warnDegenerate_(0),
     cutPoints_(0),
     cutEdges_(0),
-    facePairToVertex_(2*max(surf1.size(), surf2.size())),
     facePairToEdge_(2*max(surf1.size(), surf2.size())),
+    facePairToEdgeId_(2*max(surf1.size(), surf2.size())),
     surf1EdgeCuts_(0),
     surf2EdgeCuts_(0)
 {
@@ -1180,8 +1329,7 @@ Foam::surfaceIntersection::surfaceIntersection
         Pout<< "surfaceIntersection : Writing intersection to intEdges.obj"
             << endl;
 
-        OFstream intStream("intEdges.obj");
-        writeOBJ(cutPoints_, cutEdges_, intStream);
+        OBJstream("intEdges.obj").write(cutEdges_, cutPoints_);
 
         // Dump all cut edges to files
         Pout<< "Dumping cut edges of surface1 to surf1EdgeCuts.obj" << endl;
@@ -1193,37 +1341,8 @@ Foam::surfaceIntersection::surfaceIntersection
         writeIntersectedEdges(surf2, surf2EdgeCuts_, edge2Stream);
     }
 
-
-    // Debugging stuff
-    {
-        // Check all facePairToVertex is used.
-        labelHashSet usedPoints;
-
-        forAllConstIter(labelPairLookup, facePairToEdge_, iter)
-        {
-            const label edgeI = iter();
-            const edge& e = cutEdges_[edgeI];
-
-            usedPoints.insert(e[0]);
-            usedPoints.insert(e[1]);
-        }
-
-        forAllConstIter(labelPairLookup, facePairToVertex_, iter)
-        {
-            const label pointi = iter();
-
-            if (!usedPoints.found(pointi))
-            {
-                WarningInFunction
-                    << "Problem: cut point:" << pointi
-                    << " coord:" << cutPoints_[pointi]
-                    << " not used by any edge" << endl;
-            }
-        }
-    }
-
     // Temporaries
-    facePairToVertex_.clear();
+    facePairToEdge_.clear();
 
     // // Cleanup any duplicate cuts?
     // mergeEdges();
@@ -1244,9 +1363,9 @@ const Foam::edgeList& Foam::surfaceIntersection::cutEdges() const
 }
 
 
-const Foam::labelPairLookup& Foam::surfaceIntersection::facePairToEdge() const
+const Foam::labelPairLookup& Foam::surfaceIntersection::facePairToEdgeId() const
 {
-    return facePairToEdge_;
+    return facePairToEdgeId_;
 }
 
 
@@ -1349,9 +1468,9 @@ void Foam::surfaceIntersection::mergeEdges()
     // if (nUniqEdges < cutEdges_.size())
     // {
     //     // Additional safety, in case the edge was replaced?
-    //     forAllIter(labelPairLookup, facePairToEdge_, iter)
+    //     forAllIters(facePairToEdge_, iter)
     //     {
-    //         iter() = edgeNumbering[iter()];
+    //         iter.object() = edgeNumbering[iter.object()];
     //     }
     // }
 
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
index 5c60958675..93e15669ee 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
@@ -50,7 +50,7 @@ Description
         Property       | Description                    | Type   | Default value
         tolerance      | Edge-length tolerance          | scalar | 1e-3
         allowEdgeHits  | Edge-end cuts another edge     | bool   | true
-        avoidDuplicates | Reduce the number of duplicate points    | bool | true
+        snap           | Snap near end-points           | bool   | true
         warnDegenerate | Number of warnings about degenerate edges | label | 0
     \endtable
 
@@ -66,11 +66,11 @@ SourceFiles
 
 #include "DynamicList.H"
 #include "point.H"
-#include "edge.H"
-#include "labelPairHashes.H"
-#include "typeInfo.H"
+#include "edgeHashes.H"
 #include "edgeList.H"
+#include "labelPairHashes.H"
 #include "pointIndexHit.H"
+#include "typeInfo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -82,6 +82,14 @@ class triSurfaceSearch;
 class triSurface;
 class edgeIntersections;
 
+
+//- Key is non-commutative pair of labels. Value is commutative pair of labels
+typedef LabelPairMap<edge> labelPairEdgeLookup;
+
+//- Map from edge back to all parents (pairs of faces)
+typedef EdgeMap<labelPairHashSet> edgelabelPairHashLookup;
+
+
 /*---------------------------------------------------------------------------*\
                      Class surfaceIntersection Declaration
 \*---------------------------------------------------------------------------*/
@@ -104,8 +112,8 @@ class surfaceIntersection
         //- Allow edge-ends to cut another edge.
         bool allowEdgeHits_;
 
-        //- Avoid creating duplicate cuts near edge ends
-        bool avoidDuplicates_;
+        //- Snap cut points near edge ends (default: true)
+        bool snapToEnd_;
 
         //- Maximum number of warnings about degenerate edges
         label warnDegenerate_;
@@ -117,27 +125,29 @@ class surfaceIntersection
         //  Reference into cutPoints.
         edgeList cutEdges_;
 
-        //- From face on surf1 and face on surf2 to intersection point
-        //  (label in cutPoints)
-        labelPairLookup facePairToVertex_;
-
         //- From face on surf1 and face on surf2 to intersection edge
-        //  (label in cutEdges)
-        labelPairLookup facePairToEdge_;
+        labelPairEdgeLookup facePairToEdge_;
 
-        //- Edges on surf1 that are cut. From edge on surf1 to label in cutPoint
-        //  If multiple cuts:sorted from edge.start to edge.end
+        //- From face on surf1 and face on surf2 to intersection edgeId
+        //  (label in cutEdges)
+        labelPairLookup facePairToEdgeId_;
+
+        //- Edges on surf1 that are cut.
+        //  From edgeId on surf1 to location in cutPoint
         labelListList surf1EdgeCuts_;
 
-        //- Edges on surf2 that are cut. From edge on surf2 to label in cutPoint
-        //  If multiple cuts:sorted from edge.start to edge.end
+        //- Edges on surf2 that are cut.
+        //  From edgeId on surf2 to location in cutPoint
         labelListList surf2EdgeCuts_;
 
         //- Temporary storage to manage edge-edge self-intersections.
-        HashSet<edge, Hash<edge>> edgeEdgeIntersection_;
+        edgeHashSet edgeEdgeIntersection_;
 
         //- Temporary storage to manage cuts/intersections from the edge ends
-        Map<label> edgeEndAsCut_;
+        Map<label> snappedEnds_;
+
+        //- Temporary storage
+        EdgeMap<label> edgeToId_;
 
 
     // Private Member Functions
@@ -145,17 +155,6 @@ class surfaceIntersection
         //- Adjust intersection options according to the dictionary entries
         void setOptions(const dictionary& dict);
 
-        //- Write points in obj format
-        static void writeOBJ(const List<point>& pts, Ostream& os);
-
-        //- Write points and edges in obj format
-        static void writeOBJ
-        (
-            const List<point>& pts,
-            const List<edge>& edges,
-            Ostream& os
-        );
-
         //- Transfer contents of List<DynamicList<..>> to List<List<..>>
         template<class T>
         static void transfer(List<DynamicList<T>>&,  List<List<T>>&);
@@ -201,8 +200,8 @@ class surfaceIntersection
             const UList<point>& points
         );
 
-        //- Update reference between faceA and faceB. Updates facePairToVertex_
-        //  (first occurrence of face pair) and facePairToEdge_ (second occ.)
+        //- Update reference between faceA and faceB.
+        //  Updates facePairToEdge_ and facePairToEdgeId_ (on the second hit)
         void storeIntersection
         (
             const enum originatingType cutFrom,
@@ -241,6 +240,9 @@ class surfaceIntersection
             List<DynamicList<label>>& surfEdgeCuts
         );
 
+        //- Join disconnected intersection points
+        void joinDisconnected(DynamicList<edge>& allCutEdges);
+
 
 public:
 
@@ -292,7 +294,7 @@ public:
         const edgeList& cutEdges() const;
 
         //- Lookup of pairs of faces to created edges
-        const labelPairLookup& facePairToEdge() const;
+        const labelPairLookup& facePairToEdgeId() const;
 
         //- Access either surf1EdgeCuts (isFirstSurface = true) or
         //  surf2EdgeCuts
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
index 6e261a86cb..d8797736b0 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
@@ -29,13 +29,14 @@ License
 #include "labelPairHashes.H"
 #include "OFstream.H"
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-void Foam::surfaceIntersection::writeOBJ
-(
-    const List<point>& pts,
-    Ostream& os
-)
+namespace Foam
+{
+
+// file-scope
+// Write points in obj format
+static void writeObjPoints(const UList<point>& pts, Ostream& os)
 {
     forAll(pts, i)
     {
@@ -44,23 +45,10 @@ void Foam::surfaceIntersection::writeOBJ
     }
 }
 
+} // End namespace Foam
 
-void Foam::surfaceIntersection::writeOBJ
-(
-    const List<point>& pts,
-    const List<edge>& edges,
-    Ostream& os
-)
-{
-    writeOBJ(pts, os);
 
-    forAll(edges, i)
-    {
-        const edge& e = edges[i];
-
-        os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
-    }
-}
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 
 // Get minimum length of all edges connected to point
@@ -254,8 +242,8 @@ void Foam::surfaceIntersection::writeIntersectedEdges
     // Dump all points (surface followed by cutPoints)
     const pointField& pts = surf.localPoints();
 
-    writeOBJ(pts, os);
-    writeOBJ(cutPoints(), os);
+    writeObjPoints(pts, os);
+    writeObjPoints(cutPoints(), os);
 
     forAll(edgeCutVerts, edgeI)
     {