diff --git a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C
index 61829b404b..1ffeb97f80 100644
--- a/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C
+++ b/applications/utilities/surface/surfaceHookUp/surfaceHookUp.C
@@ -201,6 +201,71 @@ void createBoundaryEdgeTrees
 }
 
 
+class findNearestOpSubset
+{
+    const indexedOctree<treeDataEdge>& tree_;
+
+    DynamicList<label>& shapeMask_;
+
+public:
+
+    findNearestOpSubset
+    (
+        const indexedOctree<treeDataEdge>& tree,
+        DynamicList<label>& shapeMask
+    )
+    :
+        tree_(tree),
+        shapeMask_(shapeMask)
+    {}
+
+    void operator()
+    (
+        const labelUList& indices,
+        const point& sample,
+
+        scalar& nearestDistSqr,
+        label& minIndex,
+        point& nearestPoint
+    ) const
+    {
+        const treeDataEdge& shape = tree_.shapes();
+
+        forAll(indices, i)
+        {
+            const label index = indices[i];
+            const label edgeIndex = shape.edgeLabels()[index];
+
+            if
+            (
+                !shapeMask_.empty()
+             && findIndex(shapeMask_, edgeIndex) != -1
+            )
+            {
+                continue;
+            }
+
+            const edge& e = shape.edges()[edgeIndex];
+
+            pointHit nearHit = e.line(shape.points()).nearestDist(sample);
+
+            // Only register hit if closest point is not an edge point
+            if (nearHit.hit())
+            {
+                scalar distSqr = sqr(nearHit.distance());
+
+                if (distSqr < nearestDistSqr)
+                {
+                    nearestDistSqr = distSqr;
+                    minIndex = index;
+                    nearestPoint = nearHit.rawPoint();
+                }
+            }
+        }
+    }
+};
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -226,7 +291,8 @@ int main(int argc, char *argv[])
     const IOdictionary dict(dictIO);
 
     const scalar dist(args.argRead<scalar>(1));
-    const scalar matchTolerance(SMALL);
+    const scalar matchTolerance(max(1e-6*dist, SMALL));
+    const label maxIters = 100;
 
     Info<< "Hooking distance = " << dist << endl;
 
@@ -289,7 +355,7 @@ int main(int argc, char *argv[])
     }
 
     label nChanged = 0;
-    label nIters = 0;
+    label nIters = 1;
 
     do
     {
@@ -339,7 +405,7 @@ int main(int argc, char *argv[])
                         (
                             samplePt,
                             sqr(dist),
-                            treeDataEdge::findNearestOpSubset
+                            findNearestOpSubset
                             (
                                 bEdgeTree,
                                 shapeMask
@@ -439,14 +505,18 @@ int main(int argc, char *argv[])
 
                     nChanged++;
 
+                    label newPointI = -1;
+
                     // Keep the points in the same place and move the edge
-                    newPoints[hitSurfI].append(newPoints[surfI][pointI]);
-
-                    // Move the points to the edges
-                    //newPoints[pointI] = eHit.hitPoint();
-                    //newPoints.append(eHit.hitPoint());
-
-                    visitedFace[hitSurfI][faceI] = true;
+                    if (hitSurfI == surfI)
+                    {
+                        newPointI = pointI;
+                    }
+                    else
+                    {
+                        newPoints[hitSurfI].append(newPoints[surfI][pointI]);
+                        newPointI = newPoints[hitSurfI].size() - 1;
+                    }
 
                     // Split the other face.
                     greenRefine
@@ -454,22 +524,23 @@ int main(int argc, char *argv[])
                         hitSurf,
                         faceI,
                         eIndex,
-                        newPoints[hitSurfI].size() - 1,
+                        newPointI,
                         newFacesFromSplit
                     );
 
+                    visitedFace[hitSurfI][faceI] = true;
+
                     forAll(newFacesFromSplit, newFaceI)
                     {
+                        const labelledTri& fN = newFacesFromSplit[newFaceI];
+
                         if (newFaceI == 0)
                         {
-                            newFaces[hitSurfI][faceI] = newFacesFromSplit[0];
+                            newFaces[hitSurfI][faceI] = fN;
                         }
                         else
                         {
-                            newFaces[hitSurfI].append
-                            (
-                                newFacesFromSplit[newFaceI]
-                            );
+                            newFaces[hitSurfI].append(fN);
                         }
                     }
                 }
@@ -502,7 +573,7 @@ int main(int argc, char *argv[])
             );
         }
 
-    } while (nChanged > 0);
+    } while (nChanged > 0 && nIters <= maxIters);
 
     Info<< endl;
 
diff --git a/src/meshTools/indexedOctree/treeDataEdge.C b/src/meshTools/indexedOctree/treeDataEdge.C
index ffd0969e8f..4d0b9b4d81 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.C
+++ b/src/meshTools/indexedOctree/treeDataEdge.C
@@ -105,17 +105,6 @@ Foam::treeDataEdge::findNearestOp::findNearestOp
 {}
 
 
-Foam::treeDataEdge::findNearestOpSubset::findNearestOpSubset
-(
-    const indexedOctree<treeDataEdge>& tree,
-    DynamicList<label>& shapeMask
-)
-:
-    tree_(tree),
-    shapeMask_(shapeMask)
-{}
-
-
 Foam::treeDataEdge::findIntersectOp::findIntersectOp
 (
     const indexedOctree<treeDataEdge>& tree
@@ -281,105 +270,6 @@ void Foam::treeDataEdge::findNearestOp::operator()
 }
 
 
-void Foam::treeDataEdge::findNearestOpSubset::operator()
-(
-    const labelUList& indices,
-    const point& sample,
-
-    scalar& nearestDistSqr,
-    label& minIndex,
-    point& nearestPoint
-) const
-{
-    const treeDataEdge& shape = tree_.shapes();
-
-    forAll(indices, i)
-    {
-        const label index = indices[i];
-        const label edgeIndex = shape.edgeLabels()[index];
-
-        if (!shapeMask_.empty() && findIndex(shapeMask_, edgeIndex) != -1)
-        {
-            continue;
-        }
-
-        const edge& e = shape.edges()[edgeIndex];
-
-        pointHit nearHit = e.line(shape.points()).nearestDist(sample);
-
-        scalar distSqr = sqr(nearHit.distance());
-
-        if (distSqr < nearestDistSqr)
-        {
-            nearestDistSqr = distSqr;
-            minIndex = index;
-            nearestPoint = nearHit.rawPoint();
-        }
-    }
-}
-
-
-void Foam::treeDataEdge::findNearestOpSubset::operator()
-(
-    const labelUList& indices,
-    const linePointRef& ln,
-
-    treeBoundBox& tightest,
-    label& minIndex,
-    point& linePoint,
-    point& nearestPoint
-) const
-{
-    const treeDataEdge& shape = tree_.shapes();
-
-    // Best so far
-    scalar nearestDistSqr = magSqr(linePoint - nearestPoint);
-
-    forAll(indices, i)
-    {
-        const label index = indices[i];
-        const label edgeIndex = shape.edgeLabels()[index];
-
-        if (!shapeMask_.empty() && findIndex(shapeMask_, edgeIndex) != -1)
-        {
-            continue;
-        }
-
-        const edge& e = shape.edges()[edgeIndex];
-
-        // Note: could do bb test ? Worthwhile?
-
-        // Nearest point on line
-        point ePoint, lnPt;
-        scalar dist = e.line(shape.points()).nearestDist(ln, ePoint, lnPt);
-        scalar distSqr = sqr(dist);
-
-        if (distSqr < nearestDistSqr)
-        {
-            nearestDistSqr = distSqr;
-            minIndex = index;
-            linePoint = lnPt;
-            nearestPoint = ePoint;
-
-            {
-                point& minPt = tightest.min();
-                minPt = min(ln.start(), ln.end());
-                minPt.x() -= dist;
-                minPt.y() -= dist;
-                minPt.z() -= dist;
-            }
-            {
-                point& maxPt = tightest.max();
-                maxPt = max(ln.start(), ln.end());
-                maxPt.x() += dist;
-                maxPt.y() += dist;
-                maxPt.z() += dist;
-            }
-        }
-    }
-}
-
-
 bool Foam::treeDataEdge::findIntersectOp::operator()
 (
     const label index,
diff --git a/src/meshTools/indexedOctree/treeDataEdge.H b/src/meshTools/indexedOctree/treeDataEdge.H
index 5a804d5f8e..12fd0d6aa1 100644
--- a/src/meshTools/indexedOctree/treeDataEdge.H
+++ b/src/meshTools/indexedOctree/treeDataEdge.H
@@ -118,43 +118,6 @@ public:
         ) const;
     };
 
-    class findNearestOpSubset
-    {
-        const indexedOctree<treeDataEdge>& tree_;
-
-        DynamicList<label>& shapeMask_;
-
-    public:
-
-        findNearestOpSubset
-        (
-            const indexedOctree<treeDataEdge>& tree,
-            DynamicList<label>& shapeMask
-        );
-
-        void operator()
-        (
-            const labelUList& indices,
-            const point& sample,
-
-            scalar& nearestDistSqr,
-            label& minIndex,
-            point& nearestPoint
-        ) const;
-
-        //- Calculates nearest (to line) point in shape.
-        //  Returns point and distance (squared)
-        void operator()
-        (
-            const labelUList& indices,
-            const linePointRef& ln,
-
-            treeBoundBox& tightest,
-            label& minIndex,
-            point& linePoint,
-            point& nearestPoint
-        ) const;
-    };
 
     class findIntersectOp
     {