diff --git a/applications/test/surfaceIntersection/Make/files b/applications/test/surfaceIntersection/Make/files
new file mode 100644
index 0000000000..e0516a472c
--- /dev/null
+++ b/applications/test/surfaceIntersection/Make/files
@@ -0,0 +1,3 @@
+Test-surfaceIntersection.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-surfaceIntersection
diff --git a/applications/test/surfaceIntersection/Make/options b/applications/test/surfaceIntersection/Make/options
new file mode 100644
index 0000000000..6ae6c04df1
--- /dev/null
+++ b/applications/test/surfaceIntersection/Make/options
@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/edgeMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools -ledgeMesh
diff --git a/applications/test/surfaceIntersection/Test-surfaceIntersection.C b/applications/test/surfaceIntersection/Test-surfaceIntersection.C
new file mode 100644
index 0000000000..8c0f83f641
--- /dev/null
+++ b/applications/test/surfaceIntersection/Test-surfaceIntersection.C
@@ -0,0 +1,234 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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/>.
+
+Application
+    Test-surfaceIntersection
+
+Description
+    Test surface-surface intersection
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "triSurface.H"
+#include "triSurfaceMesh.H"
+#include "surfaceIntersection.H"
+#include "OFstream.H"
+
+using namespace Foam;
+
+autoPtr<triSurface> loadSurface
+(
+    const Foam::Time& runTime,
+    const fileName& surfName
+)
+{
+    Info<< "Reading surface " << surfName << endl;
+
+    const fileName fallback =
+        runTime.constantPath()/triSurfaceMesh::meshSubDir/surfName;
+
+    autoPtr<triSurface> surfPtr;
+    if (isFile(surfName))
+    {
+        surfPtr.set(new triSurface(surfName));
+    }
+    else if (isFile(fallback))
+    {
+        surfPtr.set(new triSurface(fallback));
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "No such file:" << surfName << exit(FatalError);
+    }
+
+    return surfPtr;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Intersection of two surfaces. Writes obj file"
+    );
+    argList::addBoolOption
+    (
+        "debug2",
+        "set surfaceIntersection debug=2"
+    );
+    argList::addBoolOption
+    (
+        "debug4",
+        "set surfaceIntersection debug=4"
+    );
+    argList::addBoolOption
+    (
+        "print",
+        "print information about cuts, etc"
+    );
+    argList::addBoolOption
+    (
+        "mergeEdges",
+        "merge duplicate edges"
+    );
+    argList::addOption
+    (
+        "mergePoints",
+        "mergeTol",
+        "merge points (and edges) using the specified tolerance"
+    );
+
+    #include "addDictOption.H"
+
+    argList::addNote
+    (
+        "test intersect of two surfaces. Writes obj file"
+    );
+    argList::noParallel();
+    argList::noFunctionObjects();
+    argList::validArgs.append("surface file");
+    argList::validArgs.append("surface file");
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    const word outputFile(args.executable() + ".obj");
+
+    const fileName surf1Name(args[1]);
+    triSurface surf1 = loadSurface(runTime, surf1Name)();
+    Info<< surf1Name << " statistics:" << endl;
+    surf1.writeStats(Info);
+    Info<< endl;
+
+    const fileName surf2Name(args[2]);
+    triSurface surf2 = loadSurface(runTime, surf2Name)();
+    Info<< surf2Name << " statistics:" << endl;
+    surf2.writeStats(Info);
+    Info<< endl;
+
+
+    if (args.optionFound("debug2"))
+    {
+        surfaceIntersection::debug |= 2;
+    }
+    if (args.optionFound("debug4"))
+    {
+        surfaceIntersection::debug |= 4;
+    }
+    const bool optPrint = args.optionFound("print");
+
+    dictionary intersectOptions;
+    if (args.optionFound("dict"))
+    {
+        intersectOptions = IOdictionary
+        (
+            IOobject
+            (
+                args["dict"],
+                runTime,
+                IOobject::MUST_READ_IF_MODIFIED,
+                IOobject::NO_WRITE
+            )
+        );
+    }
+
+    intersectOptions.writeEntry("intersectOptions", Info);
+    Info<< endl;
+
+    triSurfaceSearch query1(surf1);
+    triSurfaceSearch query2(surf2);
+    surfaceIntersection cuts(query1, query2, intersectOptions);
+
+    Info<<"intersection "
+        << cuts.cutPoints().size() << " points "
+        << cuts.cutEdges().size()  << " edges" << nl;
+
+    if (optPrint)
+    {
+        Info<< "surf1-cuts: " << cuts.surf1EdgeCuts() << nl
+            << "surf2-cuts: " << cuts.surf2EdgeCuts() << nl
+            << "face-pairs: " << cuts.facePairToEdge() << nl
+            << "edges: " << cuts.cutEdges() << nl;
+    }
+
+    word mergeOp;
+    if (args.optionFound("mergePoints"))
+    {
+        cuts.mergePoints(args.optionRead<scalar>("mergePoints"));
+        mergeOp = "mergePoints";
+    }
+    else if (args.optionFound("mergeEdges"))
+    {
+        cuts.mergeEdges();
+        mergeOp = "mergeEdges";
+    }
+
+    if (!mergeOp.empty())
+    {
+        Info<< mergeOp << ": "
+            << cuts.cutPoints().size() << " points "
+            << cuts.cutEdges().size()  << " edges" << nl;
+
+        if (optPrint)
+        {
+            Info<< "surf1-cuts: " << cuts.surf1EdgeCuts() << nl
+                << "surf2-cuts: " << cuts.surf2EdgeCuts() << nl
+                << "face-pairs: " << cuts.facePairToEdge() << nl
+                << "edges: " << cuts.cutEdges() << nl;
+        }
+    }
+
+    const pointField& points = cuts.cutPoints();
+    const edgeList&    edges = cuts.cutEdges();
+
+    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;
+        }
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceFeatureExtract/Allwclean b/applications/utilities/surface/surfaceFeatureExtract/Allwclean
new file mode 100755
index 0000000000..421b6686c4
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/Allwclean
@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+wclean libso extractionMethod
+wclean .
+
+#------------------------------------------------------------------------------
diff --git a/applications/utilities/surface/surfaceFeatureExtract/Allwmake b/applications/utilities/surface/surfaceFeatureExtract/Allwmake
new file mode 100755
index 0000000000..ee26d27fc3
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/Allwmake
@@ -0,0 +1,9 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Parse arguments for library compilation
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+
+(wmake libso extractionMethod && wmake)
+
+#------------------------------------------------------------------------------
diff --git a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/Make/files b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/Make/files
index 29dd2b73eb..8c9bb27660 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/Make/files
+++ b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/Make/files
@@ -1,6 +1,7 @@
 method = .
 $(method)/surfaceFeaturesExtraction.C
-$(method)/extractFromSurface.C
 $(method)/extractFromFile.C
+$(method)/extractFromNone.C
+$(method)/extractFromSurface.C
 
 LIB = $(FOAM_LIBBIN)/libsurfaceFeatureExtract
diff --git a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromFile.C b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromFile.C
index 437024fbac..d9fb9b1dde 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromFile.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromFile.C
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "extractFromFile.H"
-#include "ListOps.H"
 #include "edgeMesh.H"
 #include "addToRunTimeSelectionTable.H"
 
diff --git a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromNone.C b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromNone.C
new file mode 100644
index 0000000000..887e9a4d65
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromNone.C
@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           |  Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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 "extractFromNone.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFeaturesExtraction
+{
+    addNamedToRunTimeSelectionTable
+    (
+        method,
+        extractFromNone,
+        dictionary,
+        none
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::surfaceFeaturesExtraction::extractFromNone::extractFromNone
+(
+    const dictionary& dict
+)
+:
+    method()
+{
+    const dictionary& coeffDict = dict.subOrEmptyDict("extractFromNoneCoeffs");
+
+    coeffDict.readIfPresent("includedAngle", includedAngle_);
+    coeffDict.readIfPresent("geometricTestOnly", geometricTestOnly_);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::surfaceFeaturesExtraction::extractFromNone::~extractFromNone()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+Foam::autoPtr<Foam::surfaceFeatures>
+Foam::surfaceFeaturesExtraction::extractFromNone::features
+(
+    const triSurface& surf
+) const
+{
+    return autoPtr<surfaceFeatures>(new surfaceFeatures(surf));
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromNone.H b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromNone.H
new file mode 100644
index 0000000000..da80a3c624
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromNone.H
@@ -0,0 +1,91 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     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::surfaceFeaturesExtraction::extractFromNone
+
+Description
+    Run-time selectable surface feature extraction.
+
+SourceFiles
+    extractionMethod.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceFeaturesExtraction_extractFromNone_H
+#define surfaceFeaturesExtraction_extractFromNone_H
+
+#include "surfaceFeaturesExtraction.H"
+#include "dictionary.H"
+#include "Switch.H"
+#include "triSurface.H"
+#include "edgeIntersections.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace surfaceFeaturesExtraction
+{
+
+/*---------------------------------------------------------------------------*\
+         Class surfaceFeaturesExtraction::extractFromNone Declaration
+\*---------------------------------------------------------------------------*/
+
+class extractFromNone
+:
+    public method
+{
+
+public:
+
+    // Constructors
+
+        //- Construct from dictionary
+        extractFromNone(const dictionary& dict);
+
+
+    //- Destructor
+    virtual ~extractFromNone();
+
+
+    //- Extracted features from surface (no-op)
+    virtual autoPtr<surfaceFeatures> features
+    (
+        const triSurface& surf
+    ) const override;
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace surfaceFeaturesExtraction
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.C b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.C
index f590bb707a..6621d6c931 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.C
@@ -26,12 +26,6 @@ License
 #include "extractFromSurface.H"
 #include "addToRunTimeSelectionTable.H"
 
-#include "triSurface.H"
-#include "triSurfaceSearch.H"
-#include "scalarField.H"
-#include "edgeIntersections.H"
-
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
diff --git a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.H b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.H
index 980db47b14..e007a8d3b9 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.H
+++ b/applications/utilities/surface/surfaceFeatureExtract/extractionMethod/extractFromSurface.H
@@ -56,10 +56,6 @@ class extractFromSurface
 :
     public method
 {
-
-    bool selfIntersection_;
-
-
 public:
 
     // Constructors
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 9dcaae6b63..3ae006e69d 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -29,12 +29,17 @@ Group
 
 Description
     Extracts and writes surface features to file. All but the basic feature
-    extraction is WIP.
+    extraction is a work-in-progress.
 
-    Curvature calculation is an implementation of the algorithm from:
-
-      "Estimating Curvatures and their Derivatives on Triangle Meshes"
-      by S. Rusinkiewicz
+    The curvature calculation is an implementation of the algorithm from:
+    \verbatim
+        "Estimating Curvatures and their Derivatives on Triangle Meshes"
+        by S. Rusinkiewicz
+        3DPVT'04 Proceedings of the 3D Data Processing,
+        Visualization, and Transmission, 2nd International Symposium
+        Pages 486-493
+        http://gfx.cs.princeton.edu/pubs/_2004_ECA/curvpaper.pdf
+    \endverbatim
 
 \*---------------------------------------------------------------------------*/
 
@@ -60,6 +65,7 @@ Description
 #include "point.H"
 #include "triadField.H"
 #include "transform.H"
+#include "triSurfaceLoader.H"
 
 using namespace Foam;
 
@@ -367,23 +373,6 @@ triSurfacePointScalarField calcCurvature
 }
 
 
-bool edgesConnected(const edge& e1, const edge& e2)
-{
-    if
-    (
-        e1.start() == e2.start()
-     || e1.start() == e2.end()
-     || e1.end() == e2.start()
-     || e1.end() == e2.end()
-    )
-    {
-        return true;
-    }
-
-    return false;
-}
-
-
 scalar calcProximityOfFeaturePoints
 (
     const List<pointIndexHit>& hitList,
@@ -462,7 +451,7 @@ scalar calcProximityOfFeatureEdges
                     const edge& e2 = efem.edges()[pHit2.index()];
 
                     // Don't refine if the edges are connected to each other
-                    if (!edgesConnected(e1, e2))
+                    if (!e1.connects(e2))
                     {
                         scalar curDist =
                             mag(pHit1.hitPoint() - pHit2.hitPoint());
@@ -480,25 +469,12 @@ scalar calcProximityOfFeatureEdges
 
 void dumpBox(const treeBoundBox& bb, const fileName& fName)
 {
-    OFstream str(fName);
+    OFstream os(fName);
 
     Info<< "Dumping bounding box " << bb << " as lines to obj file "
-        << str.name() << endl;
+        << os.name() << endl;
 
-
-    pointField boxPoints(bb.points());
-
-    forAll(boxPoints, i)
-    {
-        meshTools::writeOBJ(str, boxPoints[i]);
-    }
-
-    forAll(treeBoundBox::edges, i)
-    {
-        const edge& e = treeBoundBox::edges[i];
-
-        str<< "l " << e[0]+1 <<  ' ' << e[1]+1 << nl;
-    }
+    meshTools::writeOBJ(os, bb);
 }
 
 
@@ -866,62 +842,6 @@ void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
 }
 
 
-// Read and combine all surfaces into a single one
-
-autoPtr<triSurface> loadSurfaces(Time& runTime, const wordList& surfNames)
-{
-    List<labelledTri> faces;
-    pointField points;
-    label regoff = 0; // region offset
-
-    forAll(surfNames, surfi)
-    {
-        const word& surfName = surfNames[surfi];
-
-        triSurface addsurf(runTime.constantPath()/"triSurface"/surfName);
-
-        List<labelledTri> addfaces(addsurf.xferFaces());
-        List<point> addpoints(addsurf.xferPoints());
-
-        if (surfi)
-        {
-            const label ptoff = points.size(); // point offset
-
-            forAll(addfaces, facei)
-            {
-                labelledTri& f = addfaces[facei];
-                forAll(f, fi)
-                {
-                    f[fi] += ptoff;
-                }
-                f.region() += regoff;
-            }
-
-            faces.append(addfaces);
-            points.append(addpoints);
-        }
-        else
-        {
-            faces.transfer(addfaces);
-            points.transfer(addpoints);
-        }
-
-        regoff += addsurf.patches().size();
-    }
-
-    return autoPtr<triSurface>
-    (
-        new triSurface
-        (
-            faces,
-            geometricSurfacePatchList(),
-            points,
-            true
-        )
-    );
-}
-
-
 int main(int argc, char *argv[])
 {
     argList::addNote
@@ -939,25 +859,15 @@ int main(int argc, char *argv[])
     const word dictName("surfaceFeatureExtractDict");
     #include "setSystemRunTimeDictionaryIO.H"
 
-    // Will be using triSurface, so filter according to what is supported
-    fileNameList validSurfaceFiles = readDir
-    (
-        runTime.path()/runTime.constant()/"triSurface",
-        fileName::FILE
-    );
-    inplaceSubsetList
-    (
-        validSurfaceFiles,
-        [](const fileName& f){ return triSurface::canRead(f); }
-    );
-
-    // sort and eliminate duplicates (eg, files with/without .gz)
-    inplaceUniqueSort(validSurfaceFiles);
-
     Info<< "Reading " << dictName << nl << endl;
-
     const IOdictionary dict(dictIO);
 
+    // Loader for available triSurface surface files
+    triSurfaceLoader loader(runTime);
+
+    // Where to write VTK output files
+    const fileName vtkOutputDir = runTime.constantPath()/"triSurface";
+
     forAllConstIter(dictionary, dict, iter)
     {
         if (!iter().isDict())
@@ -972,115 +882,88 @@ int main(int argc, char *argv[])
             continue;
         }
 
-        autoPtr<surfaceFeaturesExtraction::method> extractPtr =
+        autoPtr<surfaceFeaturesExtraction::method> extractor =
             surfaceFeaturesExtraction::method::New
             (
                 surfaceDict
             );
 
-        const surfaceFeaturesExtraction::method& extract = extractPtr();
-
-        // Output name, cleansed of extensions
-        const word sFeatFileName =
+        // The output name, cleansed of extensions
+        // Optional "output" entry, or the dictionary name.
+        const word outputName =
             fileName
             (
                 surfaceDict.lookupOrDefault<word>("output", iter().keyword())
             ).lessExt();
 
-        wordList surfFileNames;
-
+        // The "surfaces" entry is normally optional, but if the sub-dictionary
+        // is itself called "surfaces", then this becomes mandatory.
+        // This provides a simple means of handling both situations without an
+        // additional switch.
         if
         (
-            iter().keyword() == "surfaces"   // mandatory
-         || surfaceDict.found("surfaces")    // or optional
+            iter().keyword() == "surfaces"  // mandatory
+         || surfaceDict.found("surfaces")   // or optional
         )
         {
-            wordReList regexs(surfaceDict.lookup("surfaces"));
-
-            labelList selected = findStrings(regexs, validSurfaceFiles);
-            surfFileNames.setSize(selected.size());
-            forAll(selected, i)
-            {
-                surfFileNames[i] = validSurfaceFiles[selected[i]].name();
-            }
-
-            if (surfFileNames.empty())
-            {
-                FatalErrorInFunction
-                    << "No surfaces specified/found for entry: "
-                    << iter().keyword()
-                    << exit(FatalError);
-            }
+            loader.select(wordReList(surfaceDict.lookup("surfaces")));
         }
         else
         {
-            surfFileNames.setSize(1);
-            surfFileNames[0] = iter().keyword();
-
-            const fileName file
-            (
-                runTime.constantPath()/"triSurface"/surfFileNames[0]
-            );
-
-            if (!isFile(file))
-            {
-                FatalErrorInFunction
-                    << "No surface: " << file.name()
-                    << exit(FatalError);
-            }
+            loader.select(iter().keyword());
         }
 
-        // DebugVar(surfFileNames);
-        // DebugVar(sFeatFileName);
+        if (loader.selected().empty())
+        {
+            FatalErrorInFunction
+                << "No surfaces specified/found for entry: "
+                << iter().keyword() << exit(FatalError);
+        }
+        // DebugVar(loader.available());
+        // DebugVar(outputName);
+
 
         Info<< "Surfaces           : ";
-        if (surfFileNames.size() == 1)
+        if (loader.selected().size() == 1)
         {
-            Info<< surfFileNames[0] << nl;
+            Info<< loader.selected()[0] << nl;
         }
         else
         {
-            Info<< flatOutput(surfFileNames) << nl;
+            Info<< flatOutput(loader.selected()) << nl;
+        }
+        Info<< "Output             : " << outputName << nl;
+
+        // Load a single file, or load and combine multiple selected files
+        autoPtr<triSurface> surfPtr = loader.load();
+        if (!surfPtr.valid() || surfPtr().empty())
+        {
+            FatalErrorInFunction
+                << "Problem loading surface(s) for entry: "
+                << iter().keyword() << exit(FatalError);
         }
 
-        Info<< "Output             : " << sFeatFileName << nl;
+        triSurface surf = surfPtr();
 
         const Switch writeVTK = surfaceDict.lookupOrDefault<Switch>
         (
-            "writeVTK", Switch::OFF
+            "writeVTK",
+            Switch::OFF
         );
         const Switch writeObj = surfaceDict.lookupOrDefault<Switch>
         (
-            "writeObj", Switch::OFF
-        );
-
-        const Switch curvature = surfaceDict.lookupOrDefault<Switch>
-        (
-            "curvature", Switch::OFF
-        );
-        const Switch featureProximity = surfaceDict.lookupOrDefault<Switch>
-        (
-            "featureProximity", Switch::OFF
-        );
-        const Switch closeness = surfaceDict.lookupOrDefault<Switch>
-        (
-            "closeness", Switch::OFF
+            "writeObj",
+            Switch::OFF
         );
 
         Info<< "write VTK: " <<  writeVTK << nl;
 
-        Info<< nl << "Feature line extraction is only valid on closed manifold "
-            << "surfaces." << endl;
+        Info<< "Feature line extraction is only valid on closed manifold "
+            << "surfaces." << nl;
 
-
-        // Read and combine all surfaces into a single one
-
-        autoPtr<triSurface> surfPtr = loadSurfaces(runTime, surfFileNames);
-        triSurface surf = surfPtr();
-
-        Info<< "Statistics:" << endl;
+        Info<< nl << "Statistics:" << nl;
         surf.writeStats(Info);
-        Info<< endl;
+        Info<< nl;
 
         // need plain faces if outputting VTK format
         faceList faces;
@@ -1097,19 +980,19 @@ int main(int argc, char *argv[])
         // Either construct features from surface & featureAngle or read set.
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-        autoPtr<surfaceFeatures> set = extract.features(surf);
+        autoPtr<surfaceFeatures> features = extractor().features(surf);
 
         // Trim set
         // ~~~~~~~~
 
         if (surfaceDict.isDict("trimFeatures"))
         {
-            dictionary trimDict = surfaceDict.subDict("trimFeatures");
-            const scalar minLen =
-                trimDict.lookupOrAddDefault<scalar>("minLen", -GREAT);
+            const dictionary& trimDict = surfaceDict.subDict("trimFeatures");
 
+            const scalar minLen =
+                trimDict.lookupOrDefault<scalar>("minLen", 0);
             const label minElem =
-                trimDict.lookupOrAddDefault<label>("minElem", 0);
+                trimDict.lookupOrDefault<label>("minElem", 0);
 
             // Trim away small groups of features
             if (minLen > 0 || minElem > 0)
@@ -1125,7 +1008,10 @@ int main(int argc, char *argv[])
                         << minElem << endl;
                 }
 
-                set().trimFeatures(minLen, minElem, extract.includedAngle());
+                features().trimFeatures
+                (
+                    minLen, minElem, extractor().includedAngle()
+                );
             }
         }
 
@@ -1134,8 +1020,9 @@ int main(int argc, char *argv[])
         // ~~~~~~
 
         // Convert to marked edges, points
-        List<surfaceFeatures::edgeStatus> edgeStat(set().toStatus());
+        List<surfaceFeatures::edgeStatus> edgeStat(features().toStatus());
 
+        // Option: "subsetFeatures" (dictionary)
         if (surfaceDict.isDict("subsetFeatures"))
         {
             const dictionary& subsetDict = surfaceDict.subDict
@@ -1143,6 +1030,7 @@ int main(int argc, char *argv[])
                 "subsetFeatures"
             );
 
+            // Suboption: "insideBox"
             if (subsetDict.found("insideBox"))
             {
                 treeBoundBox bb(subsetDict.lookup("insideBox")());
@@ -1152,6 +1040,7 @@ int main(int argc, char *argv[])
 
                 deleteBox(surf, bb, false, edgeStat);
             }
+            // Suboption: "outsideBox"
             else if (subsetDict.found("outsideBox"))
             {
                 treeBoundBox bb(subsetDict.lookup("outsideBox")());
@@ -1186,17 +1075,15 @@ int main(int argc, char *argv[])
                         (
                             surf,
                             1e-5,   //tol,
-                            extract.includedAngle(),
+                            extractor().includedAngle(),
                             edgeI
                         );
                     }
                 }
             }
 
-            const Switch openEdges =
-                subsetDict.lookupOrDefault<Switch>("openEdges", "yes");
-
-            if (!openEdges)
+            // Suboption: "openEdges" (false: remove open edges)
+            if (!subsetDict.lookupOrDefault<bool>("openEdges", true))
             {
                 Info<< "Removing all open edges"
                     << " (edges with 1 connected face)" << endl;
@@ -1210,6 +1097,7 @@ int main(int argc, char *argv[])
                 }
             }
 
+            // Suboption: "plane"
             if (subsetDict.found("plane"))
             {
                 plane cutPlane(subsetDict.lookup("plane")());
@@ -1225,7 +1113,7 @@ int main(int argc, char *argv[])
 
 
         surfaceFeatures newSet(surf);
-        newSet.setFromStatus(edgeStat, extract.includedAngle());
+        newSet.setFromStatus(edgeStat, extractor().includedAngle());
 
         Info<< nl
             << "Initial feature set:" << nl
@@ -1263,7 +1151,7 @@ int main(int argc, char *argv[])
         (
             newSet,
             runTime,
-            sFeatFileName + ".extendedFeatureEdgeMesh",
+            outputName + ".extendedFeatureEdgeMesh",
             surfBaffleRegions
         );
 
@@ -1292,40 +1180,41 @@ int main(int argc, char *argv[])
             feMesh.add(addFeMesh);
         }
 
-
         if (surfaceDict.lookupOrDefault<bool>("selfIntersection", false))
         {
-            // TODO: perturb tolerance
+            // TODO: perturbance tolerance?
 
             triSurfaceSearch query(surf);
-            surfaceIntersection intersect
+            surfaceIntersection intersect(query, surfaceDict);
+
+            intersect.mergePoints(5*SMALL);
+
+            labelPair sizeInfo
             (
-                query,
-                surfaceDict.lookupOrDefault<scalar>
-                (
-                    "planarTolerance",
-                    surfaceIntersection::defaultTolerance
-                )
+                intersect.cutPoints().size(),
+                intersect.cutEdges().size()
             );
 
-            // surf.write("selfIntersection-input.obj");
-
-            Info<<"self-intersection "
-                << intersect.cutEdges().size() << " edges "
-                << intersect.cutPoints().size() << " points" << nl;
-
             if (intersect.cutEdges().size())
             {
                 extendedEdgeMesh addMesh
                 (
-                    xferCopy<pointField>(intersect.cutPoints()),
-                    xferCopy<edgeList>(intersect.cutEdges())
+                    intersect.cutPoints(),
+                    intersect.cutEdges()
                 );
 
+                addMesh.mergePoints(5*SMALL);
                 feMesh.add(addMesh);
+
+                sizeInfo[0] = addMesh.points().size();
+                sizeInfo[1] = addMesh.edges().size();
             }
+            Info<< "Self intersection:" << nl
+                    << "    points : " << sizeInfo[0] << nl
+                    << "    edges  : " << sizeInfo[1] << nl;
         }
 
+
         Info<< nl
             << "Final feature set:" << nl;
         writeStats(feMesh, Info);
@@ -1337,111 +1226,47 @@ int main(int argc, char *argv[])
 
         if (writeObj)
         {
-            feMesh.writeObj(feMesh.path()/sFeatFileName);
+            feMesh.writeObj(feMesh.path()/outputName);
         }
 
         feMesh.write();
 
         // Write a featureEdgeMesh for backwards compatibility
-        featureEdgeMesh bfeMesh
-        (
-            IOobject
+        if (true)
+        {
+            featureEdgeMesh bfeMesh
             (
-                sFeatFileName + ".eMesh",   // name
-                runTime.constant(),         // instance
-                "triSurface",
-                runTime,                    // registry
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE,
-                false
-            ),
-            feMesh.points(),
-            feMesh.edges()
-        );
+                IOobject
+                (
+                    outputName + ".eMesh",      // name
+                    runTime.constant(),         // instance
+                    "triSurface",
+                    runTime,                    // registry
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE,
+                    false
+                ),
+                feMesh.points(),
+                feMesh.edges()
+            );
 
-        Info<< nl << "Writing featureEdgeMesh to "
-            << bfeMesh.objectPath() << endl;
+            Info<< nl << "Writing featureEdgeMesh to "
+                << bfeMesh.objectPath() << endl;
 
-        bfeMesh.regIOobject::write();
+            bfeMesh.regIOobject::write();
+        }
 
-    // Find close features
-
-    // // Dummy trim operation to mark features
-    // labelList featureEdgeIndexing = newSet.trimFeatures(-GREAT, 0);
-
-    // scalarField surfacePtFeatureIndex(surf.points().size(), -1);
-
-    // forAll(newSet.featureEdges(), eI)
-    // {
-    //     const edge& e = surf.edges()[newSet.featureEdges()[eI]];
-
-    //     surfacePtFeatureIndex[surf.meshPoints()[e.start()]] =
-    //     featureEdgeIndexing[newSet.featureEdges()[eI]];
-
-    //     surfacePtFeatureIndex[surf.meshPoints()[e.end()]] =
-    //     featureEdgeIndexing[newSet.featureEdges()[eI]];
-    // }
-
-    // if (writeVTK)
-    // {
-    //     vtkSurfaceWriter().write
-    //     (
-    //         runTime.constant()/"triSurface",    // outputDir
-    //         sFeatFileName,                      // surfaceName
-    //         surf.points(),
-    //         faces,
-    //         "surfacePtFeatureIndex",            // fieldName
-    //         surfacePtFeatureIndex,
-    //         true,                               // isNodeValues
-    //         true                                // verbose
-    //     );
-    // }
-
-    // Random rndGen(343267);
-
-    // treeBoundBox surfBB
-    // (
-    //     treeBoundBox(searchSurf.bounds()).extend(rndGen, 1e-4)
-    // );
-
-    // surfBB.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-    // surfBB.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
-
-    // indexedOctree<treeDataEdge> ftEdTree
-    // (
-    //     treeDataEdge
-    //     (
-    //         false,
-    //         surf.edges(),
-    //         surf.localPoints(),
-    //         newSet.featureEdges()
-    //     ),
-    //     surfBB,
-    //     8,      // maxLevel
-    //     10,     // leafsize
-    //     3.0     // duplicity
-    // );
-
-    // labelList nearPoints = ftEdTree.findBox
-    // (
-    //     treeBoundBox
-    //     (
-    //         sPt - featureSearchSpan*Foam::vector::one,
-    //         sPt + featureSearchSpan*Foam::vector::one
-    //     )
-    // );
-
-        if (closeness)
+        // Option: "closeness"
+        if (surfaceDict.lookupOrDefault<bool>("closeness", false))
         {
             Info<< nl << "Extracting internal and external closeness of "
                 << "surface." << endl;
 
-
             triSurfaceMesh searchSurf
             (
                 IOobject
                 (
-                    sFeatFileName + ".closeness",
+                    outputName + ".closeness",
                     runTime.constant(),
                     "triSurface",
                     runTime,
@@ -1481,8 +1306,8 @@ int main(int argc, char *argv[])
 
             // Info<< "span " << span << endl;
 
-            pointField start(searchSurf.faceCentres() - span*normals);
-            pointField end(searchSurf.faceCentres() + span*normals);
+            const pointField start(searchSurf.faceCentres() - span*normals);
+            const pointField end(searchSurf.faceCentres() + span*normals);
             const pointField& faceCentres = searchSurf.faceCentres();
 
             List<List<pointIndexHit>> allHitInfo;
@@ -1649,7 +1474,7 @@ int main(int argc, char *argv[])
             (
                 IOobject
                 (
-                    sFeatFileName + ".internalCloseness",
+                    outputName + ".internalCloseness",
                     runTime.constant(),
                     "triSurface",
                     runTime,
@@ -1667,7 +1492,7 @@ int main(int argc, char *argv[])
             (
                 IOobject
                 (
-                    sFeatFileName + ".externalCloseness",
+                    outputName + ".externalCloseness",
                     runTime.constant(),
                     "triSurface",
                     runTime,
@@ -1685,8 +1510,8 @@ int main(int argc, char *argv[])
             {
                 vtkSurfaceWriter().write
                 (
-                    runTime.constantPath()/"triSurface",// outputDir
-                    sFeatFileName,                      // surfaceName
+                    vtkOutputDir,
+                    outputName,
                     meshedSurfRef
                     (
                         surf.points(),
@@ -1700,8 +1525,8 @@ int main(int argc, char *argv[])
 
                 vtkSurfaceWriter().write
                 (
-                    runTime.constantPath()/"triSurface",// outputDir
-                    sFeatFileName,                      // surfaceName
+                    vtkOutputDir,
+                    outputName,
                     meshedSurfRef
                     (
                         surf.points(),
@@ -1715,8 +1540,8 @@ int main(int argc, char *argv[])
             }
         }
 
-
-        if (curvature)
+        // Option: "curvature"
+        if (surfaceDict.lookupOrDefault<bool>("curvature", false))
         {
             Info<< nl << "Extracting curvature of surface at the points."
                 << endl;
@@ -1726,7 +1551,7 @@ int main(int argc, char *argv[])
 
             triSurfacePointScalarField k = calcCurvature
             (
-                sFeatFileName,
+                outputName,
                 runTime,
                 surf,
                 pointNormals,
@@ -1739,8 +1564,8 @@ int main(int argc, char *argv[])
             {
                 vtkSurfaceWriter().write
                 (
-                    runTime.constantPath()/"triSurface",// outputDir
-                    sFeatFileName,                      // surfaceName
+                    vtkOutputDir,
+                    outputName,
                     meshedSurfRef
                     (
                         surf.points(),
@@ -1754,8 +1579,8 @@ int main(int argc, char *argv[])
             }
         }
 
-
-        if (featureProximity)
+        // Option: "featureProximity"
+        if (surfaceDict.lookupOrDefault<bool>("featureProximity", false))
         {
             Info<< nl << "Extracting proximity of close feature points and "
                 << "edges to the surface" << endl;
@@ -1802,7 +1627,7 @@ int main(int argc, char *argv[])
             (
                 IOobject
                 (
-                    sFeatFileName + ".featureProximity",
+                    outputName + ".featureProximity",
                     runTime.constant(),
                     "triSurface",
                     runTime,
@@ -1820,8 +1645,8 @@ int main(int argc, char *argv[])
             {
                 vtkSurfaceWriter().write
                 (
-                    runTime.constantPath()/"triSurface",// outputDir
-                    sFeatFileName,                      // surfaceName
+                    vtkOutputDir,
+                    outputName,
                     meshedSurfRef
                     (
                         surf.points(),
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
index f110e06318..d04bbc7958 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractDict
@@ -16,7 +16,7 @@ FoamFile
 
 surface1.stl
 {
-    // How to obtain raw features (extractFromFile || extractFromSurface)
+    // How to obtain raw features (extractFromFile | extractFromSurface | none)
     extractionMethod    extractFromSurface;
 
     extractFromSurfaceCoeffs
@@ -28,19 +28,25 @@ surface1.stl
         includedAngle   120;
 
         // Do not mark region edges
-        geometricTestOnly       yes;
+        geometricTestOnly  yes;
     }
 
+    // Generate additional features from self-intersect
+    selfIntersection    false;
+
+    // Tolerance for surface intersections
+    tolerance           1e-3;
+
     // Write options
 
         // Write features to obj format for postprocessing
-        writeObj                yes;
+        writeObj        yes;
 }
 
 
 surface2.nas
 {
-    // How to obtain raw features (extractFromFile || extractFromSurface)
+    // How to obtain raw features (extractFromFile | extractFromSurface | none)
     extractionMethod    extractFromFile;
 
     extractFromFileCoeffs
@@ -91,36 +97,42 @@ surface2.nas
     }
 
     // Output the curvature of the surface
-    curvature               no;
+    curvature           no;
 
     // Output the proximity of feature points and edges to each other
-    featureProximity        no;
+    featureProximity    no;
 
     // The maximum search distance to use when looking for other feature
     // points and edges
-    maxFeatureProximity     1;
+    maxFeatureProximity 1;
 
     // Out put the closeness of surface elements to other surface elements.
-    closeness               no;
+    closeness           no;
+
+    // Generate additional features from self-intersect
+    selfIntersection    false;
+
+    // Tolerance for surface intersections
+    tolerance           1e-3;
 
     // Write options
 
         // Write features to obj format for postprocessing
-        writeObj                yes;
+        writeObj        yes;
 
         // Write surface proximity and curvature fields to vtk format
         // for postprocessing
-        writeVTK                no;
+        writeVTK        no;
 }
 
 
-// Handle multiple surfaces
+// Handle single or multiple surfaces
 //
 // - If the dictionary is named 'surfaces', it must also contain a 'surfaces'
 //   entry (wordRe list).
 //
-// - If other dictionaries may contain a 'surfaces' entry, it will be taken
-//   for the input.
+// - If other dictionaries contain a 'surfaces' entry,
+//   it will be taken for the input.
 //
 surfaces
 {
@@ -131,24 +143,30 @@ surfaces
     // Base output name (optional)
     // output              surfaces;
 
-    // Generate features from self-intersect
+    // Generate additional features from self-intersect
     selfIntersection    true;
 
-    // Tolerance for self-intersect
-    planarTolerance     1e-3;
+    // Tolerance for surface intersections
+    tolerance           1e-3;
 
     extractFromSurfaceCoeffs
     {
         includedAngle   120;
 
         // Do not mark region edges
-        geometricTestOnly       yes;
+        geometricTestOnly  yes;
+    }
+
+    extractFromNoneCoeffs
+    {
+        includedAngle   120;
+
     }
 
     // Write options
 
         // Write features to obj format for postprocessing
-        writeObj                yes;
+        writeObj        yes;
 }
 
 // ************************************************************************* //
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
index 4eca27e13d..737387a5e1 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,175 +40,32 @@ namespace Foam
 defineTypeNameAndDebug(surfaceIntersection, 0);
 }
 
-const Foam::scalar Foam::surfaceIntersection::defaultTolerance = 1e-3;
-
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Checks if there exists a special topological situation that causes
-// edge and the face it hit not to be recognized.
-//
-// For now if the face shares a point with the edge
-bool Foam::surfaceIntersection::excludeEdgeHit
-(
-    const triSurface& surf,
-    const label edgeI,
-    const label facei,
-    const scalar
-)
+void Foam::surfaceIntersection::setOptions(const dictionary& dict)
 {
-    const triSurface::FaceType& f = surf.localFaces()[facei];
-    const edge& e = surf.edges()[edgeI];
-
-    forAll(f, fp)
-    {
-        if (f[fp] == e.start() || f[fp] == e.end())
-        {
-            return true;
-        }
-    }
-
-// {
-//        // Get edge vector
-//        vector eVec = e.vec(surf.localPoints());
-//        eVec /= mag(eVec) + VSMALL;
-//
-//        const labelList& eLabels = surf.faceEdges()[facei];
-//
-//        // Get edge vector of 0th edge of face
-//        vector e0Vec = surf.edges()[eLabels[0]].vec(surf.localPoints());
-//        e0Vec /= mag(e0Vec) + VSMALL;
-//
-//        vector n = e0Vec ^ eVec;
-//
-//        if (mag(n) < SMALL)
-//        {
-//            // e0 is aligned with e. Choose next edge of face.
-//            vector e1Vec = surf.edges()[eLabels[1]].vec(surf.localPoints());
-//            e1Vec /= mag(e1Vec) + VSMALL;
-//
-//            n = e1Vec ^ eVec;
-//
-//            if (mag(n) < SMALL)
-//            {
-//                // Problematic triangle. Two edges aligned with edgeI. Give
-//                // up.
-//                return true;
-//            }
-//        }
-//
-//        // Check if same as faceNormal
-//        if (mag(n & surf.faceNormals()[facei]) > 1-tol)
-//        {
-//
-//            Pout<< "edge:" << e << "  face:" << facei
-//                << "  e0Vec:" << e0Vec << "  n:" << n
-//                << "  normalComponent:" << (n & surf.faceNormals()[facei])
-//                << "  tol:" << tol << endl;
-//
-//            return true;
-//        }
-//        else
-//        {
-//            return false;
-//        }
-// }
-
-    return false;
+    dict.readIfPresent("tolerance",       tolerance_);
+    dict.readIfPresent("allowEdgeHits",   allowEdgeHits_);
+    dict.readIfPresent("avoidDuplicates", avoidDuplicates_);
+    dict.readIfPresent("warnDegenerate",  warnDegenerate_);
 }
 
 
-//// Find intersection of plane with edges of hitFacei. Returns
-//// - edgeI
-//// - intersection point
-//Foam::pointIndexHit Foam::surfaceIntersection::faceEdgeIntersection
-//(
-//    const triSurface& surf,
-//    const label hitFacei,
-//
-//    const vector& n,
-//    const point& eStart,
-//    const point& eEnd
-//)
-//{
-//    pointIndexHit pInter;
-//
-//    const pointField& points = surf.points();
-//
-//    const triSurface::FaceType& f = surf.localFaces()[hitFacei];
-//
-//    // Plane for intersect test.
-//    plane pl(eStart, n);
-//
-//    forAll(f, fp)
-//    {
-//        label fp1 = f.fcIndex(fp);
-//
-//        const point& start = points[f[fp]];
-//        const point& end = points[f[fp1]];
-//
-//        vector eVec(end - start);
-//
-//        scalar s = pl.normalIntersect(start, eVec);
-//
-//        if (s < 0 || s > 1)
-//        {
-//            pInter.setPoint(start + s*eVec);
-//
-//            // Check if is correct one: orientation walking
-//            //  eStart - eEnd - hitPoint should be opposite n
-//            vector n2(triPointRef(start, end, pInter.hitPoint()).normal());
-//
-//            Pout<< "plane normal:" << n
-//                << "  start:" << start << "  end:" << end
-//                << "  hit at:" << pInter.hitPoint()
-//                << "  resulting normal:" << n2 << endl;
-//
-//            if ((n2 & n) < 0)
-//            {
-//                pInter.setHit();
-//
-//                // Find corresponding edge between f[fp] f[fp1]
-//                label edgeI =
-//                    meshTools::findEdge
-//                    (
-//                        surf.edges(),
-//                        surf.faceEdges()[hitFacei],
-//                        f[fp],
-//                        f[fp1]
-//                    );
-//
-//                pInter.setIndex(edgeI);
-//
-//                return pInter;
-//            }
-//        }
-//    }
-//
-//    FatalErrorInFunction
-//        << "Did not find intersection of plane " << pl
-//        << " with edges of face " << hitFacei << " verts:" << f
-//        << abort(FatalError);
-//
-//    return pInter;
-//}
-
-
 void Foam::surfaceIntersection::storeIntersection
 (
     const enum originatingType cutFrom,
     const labelList& facesA,
     const label faceB,
-    DynamicList<edge>& allCutEdges,
-    DynamicList<point>& allCutPoints
+    const UList<point>& allCutPoints,
+    const label cutPointId,
+    DynamicList<edge>& allCutEdges
 )
 {
     // Our lookup for two faces - populate with faceB (invariant)
     // Normally always have face from the first surface as first element
     labelPair twoFaces(faceB, faceB);
 
-    const label pointId = allCutPoints.size()-1;
-
     forAll(facesA, facesAI)
     {
         const label faceA = facesA[facesAI];
@@ -249,34 +106,82 @@ void Foam::surfaceIntersection::storeIntersection
         if (iter == facePairToVertex_.end())
         {
             // New intersection. Store face-face intersection.
-            facePairToVertex_.insert(twoFaces, pointId);
+            if (debug & 4)
+            {
+                Pout<< "intersect faces " << twoFaces
+                    << " point-1: " << cutPointId << " = "
+                    << allCutPoints[cutPointId] << endl;
+            }
+
+            facePairToVertex_.insert(twoFaces, cutPointId);
+        }
+        else if (*iter == cutPointId)
+        {
+            // Avoid creating an edge if cutPointId had already been used
+
+            if (debug & 4)
+            {
+                Pout<< "intersect faces " << twoFaces
+                    << " dup-point: " << cutPointId << endl;
+            }
         }
         else
         {
+            const label nextEdgeId = allCutEdges.size();
+            const edge nextEdge(*iter, cutPointId, true);
+
             // 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
-            const point& prevHit = allCutPoints[*iter];
-            const point& thisHit = allCutPoints.last();
-
-            if (mag(prevHit - thisHit) < SMALL)
+            if (nextEdge.mag(allCutPoints) < SMALL)
             {
-                WarningInFunction
-                    << "Encountered degenerate edge between face "
-                    << twoFaces[0] << " on first surface and face "
-                    << twoFaces[1] << " on second surface" << endl
-                    << "Point on first surface:" << prevHit << endl
-                    << "Point on second surface:" << thisHit << endl
-                    << endl;
+                // Don't normally emit warnings, since these also arise for
+                // manifold connections. For example,
+                //
+                //   e1|  /e2
+                //     | /
+                //     |/
+                // ----.---- plane
+                //
+                // 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)
+                {
+                    --warnDegenerate_;
+                    WarningInFunction
+                        << "Degenerate edge between faces " << twoFaces
+                        << " on 1st/2nd surface with points "
+                        << nextEdge.line(allCutPoints)
+                        << endl;
+                }
+                else if (debug & 4)
+                {
+                    Pout<< "degenerate edge face-pair " << twoFaces << " "
+                        << *iter << " point " << allCutPoints[*iter]
+                        << endl;
+                }
             }
-            else
+            else if (facePairToEdge_.insert(twoFaces, nextEdgeId))
             {
-                allCutEdges.append(edge(*iter, pointId));
+                // Record complete (line) intersection of two faces
 
-                // Record intersection of faces on surf
-                facePairToEdge_.insert(twoFaces, allCutEdges.size()-1);
+                allCutEdges.append(nextEdge);
+
+                if (debug & 4)
+                {
+                    Pout<< "create edge - faces " << twoFaces << " edge#"
+                        << nextEdgeId << " edge " << nextEdge
+                        << " = " << nextEdge.line(allCutPoints)
+                        << endl;
+                }
             }
         }
     }
@@ -302,8 +207,8 @@ void Foam::surfaceIntersection::classifyHit
     const label edgeI,
     const pointIndexHit& pHit,
 
-    DynamicList<edge>& allCutEdges,
     DynamicList<point>& allCutPoints,
+    DynamicList<edge>& allCutEdges,
     List<DynamicList<label>>& surfEdgeCuts
 )
 {
@@ -317,6 +222,7 @@ void Foam::surfaceIntersection::classifyHit
     // Classify point on surface2
 
     const triSurface::FaceType& f2 = surf2.localFaces()[surf2Facei];
+    const pointField& surf1Pts = surf1.localPoints();
     const pointField& surf2Pts = surf2.localPoints();
 
     label nearType, nearLabel;
@@ -331,7 +237,7 @@ void Foam::surfaceIntersection::classifyHit
             surf1PointTol[e.end()],
             pHit.hitPoint(),
             e,
-            surf1.localPoints()
+            surf1Pts
         );
 
     if (nearType == triPointRef::POINT)
@@ -341,28 +247,71 @@ void Foam::surfaceIntersection::classifyHit
             // 1. Point hits point. Do nothing.
             if (debug & 2)
             {
-                Pout<< pHit.hitPoint() << " is surf1:"
-                    << " end point of edge " << e
+                Pout<< "hit-type[1] " << pHit.hitPoint() << " is surf1:"
+                    << " end point of edge[" << edgeI << "] " << e
+                    << "==" << e.line(surf1Pts)
                     << " surf2: vertex " << f2[nearLabel]
-                    << " coord:" << surf2Pts[f2[nearLabel]] << endl;
+                    << " coord:" << surf2Pts[f2[nearLabel]]
+                    << " - suppressed" << endl;
             }
         }
         else
         {
             // 2. Edge hits point. Cut edge with new point.
-            if (debug & 2)
+            bool cached = false;
+            label cutPointId = -1;
+            const label nearVert = f2[nearLabel];
+
+            // For self-intersection, we have tolerances for each point
+            // (surf2 is actually surf1) so we shift the hit to coincide
+            // identically.
+            if (cutFrom == surfaceIntersection::SELF)
             {
-                Pout<< pHit.hitPoint() << " is surf1:"
-                    << " somewhere on edge " << e
-                    << " surf2: vertex " << f2[nearLabel]
-                    << " coord:" << surf2Pts[f2[nearLabel]] << endl;
+                const point& nearPt = surf1Pts[nearVert];
+
+                if (mag(pHit.hitPoint() - nearPt) < surf1PointTol[nearVert])
+                {
+                    cutPointId = allCutPoints.size();
+
+                    if (avoidDuplicates_)
+                    {
+                        if (edgeEndAsCut_.insert(nearVert, cutPointId))
+                        {
+                            // First time with this end-point
+                            allCutPoints.append(nearPt);
+                        }
+                        else
+                        {
+                            // Already seen this end point
+                            cutPointId = edgeEndAsCut_[nearVert];
+                            cached = true;
+                        }
+                    }
+                    else
+                    {
+                        allCutPoints.append(nearPt);
+                    }
+                }
             }
 
-            allCutPoints.append(pHit.hitPoint());
-            surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
+            if (debug & 2)
+            {
+                Pout<< "hit-type[2] " << pHit.hitPoint() << " is surf1:"
+                    << " from edge[" << edgeI << "] " << e
+                    << " surf2: vertex " << f2[nearLabel]
+                    << " coord:" << surf2Pts[f2[nearLabel]]
+                    << " - "
+                    << (cached ? "cached" : "stored") << endl;
+            }
+
+            if (cutPointId == -1)
+            {
+                cutPointId = allCutPoints.size();
+                allCutPoints.append(pHit.hitPoint());
+            }
+            surfEdgeCuts[edgeI].append(cutPointId);
 
             const labelList& facesB = surf2.pointFaces()[f2[nearLabel]];
-
             forAll(facesB, faceBI)
             {
                 storeIntersection
@@ -370,8 +319,9 @@ void Foam::surfaceIntersection::classifyHit
                     cutFrom,
                     facesA,
                     facesB[faceBI],
-                    allCutEdges,
-                    allCutPoints
+                    allCutPoints,
+                    cutPointId,
+                    allCutEdges
                 );
             }
         }
@@ -380,18 +330,114 @@ void Foam::surfaceIntersection::classifyHit
     {
         if (edgeEnd >= 0)
         {
-            // 3. Point hits edge. Do nothing on this side. Reverse
-            // is handled by 2 (edge hits point)
-            label edge2I = getEdge(surf2, surf2Facei, nearLabel);
+            // 3. Point hits edge.
+            // Normally do nothing on this side since the reverse
+            // (edge hits point) is handled by 2.
+            // However, if the surfaces are separated by a minor gap,
+            // the end-point of a tolerance-extended edge can intersect another
+            // edge without itself being intersected by an edge.
+
+            const label edge2I = getEdge(surf2, surf2Facei, nearLabel);
             const edge& e2 = surf2.edges()[edge2I];
+            const label nearVert  = (edgeEnd == 0 ? e.start() : e.end());
+
+            label cutPointId = -1;
+
+            // Storage treatment
+            // =0: nothing/ignore
+            // >0: store point/edge-cut. Attempt to create new edge.
+            // <0: store point/edge-cut only
+            int handling = (allowEdgeHits_ ? 1 : 0);
+            if (allowEdgeHits_ && cutFrom == surfaceIntersection::SELF)
+            {
+                // The edge-edge intersection is hashed as an 'edge' to
+                // exploit the commutative lookup.
+                // Ie, only do the cut once
+                const edge intersect(edgeI, edge2I);
+
+                if (e2.found(nearVert))
+                {
+                    // Actually the same as #1 above, but missed due to
+                    // tolerancing
+                    handling = 0; // suppress
+                }
+                else if (edgeEdgeIntersection_.insert(intersect))
+                {
+                    const point& nearPt = surf1Pts[nearVert];
+
+                    if
+                    (
+                        mag(pHit.hitPoint() - nearPt) < surf1PointTol[nearVert]
+                    )
+                    {
+                        cutPointId = allCutPoints.size();
+
+                        if (avoidDuplicates_)
+                        {
+                            if (edgeEndAsCut_.insert(nearVert, cutPointId))
+                            {
+                                // First time with this end-point
+                                allCutPoints.append(nearPt);
+                            }
+                            else
+                            {
+                                // Already seen this end point
+                                cutPointId = edgeEndAsCut_[nearVert];
+                                handling = 2;  // cached
+                            }
+                        }
+                        else
+                        {
+                            allCutPoints.append(nearPt);
+                        }
+                    }
+                }
+                else
+                {
+                    handling = 0; // ignore - already did this interaction
+                }
+            }
 
             if (debug & 2)
             {
-                Pout<< pHit.hitPoint() << " is surf1:"
-                    << " end point of edge " << e
-                    << " surf2: edge " << e2
-                    << " coords:" << surf2Pts[e2.start()]
-                    << surf2Pts[e2.end()] << endl;
+                Pout<< "hit-type[3] " << pHit.hitPoint() << " is surf1:"
+                    << " end point of edge[" << edgeI << "] " << e
+                    << "==" << e.line(surf1Pts)
+                    << " surf2: edge[" << edge2I << "] " << e2
+                    << " coords:" << e2.line(surf2Pts)
+                    << " - "
+                    << (
+                           handling > 1
+                         ? "cached" : handling
+                         ? "stored" : "suppressed"
+                       ) << endl;
+            }
+
+            if (handling)
+            {
+                if (cutPointId == -1)
+                {
+                    cutPointId = allCutPoints.size();
+                    allCutPoints.append(pHit.hitPoint());
+                }
+                surfEdgeCuts[edgeI].append(cutPointId);
+            }
+
+            if (handling > 0)
+            {
+                const labelList& facesB = surf2.edgeFaces()[edge2I];
+                forAll(facesB, faceBI)
+                {
+                    storeIntersection
+                    (
+                        cutFrom,
+                        facesA,
+                        facesB[faceBI],
+                        allCutPoints,
+                        cutPointId,
+                        allCutEdges
+                    );
+                }
             }
         }
         else
@@ -402,32 +448,108 @@ void Foam::surfaceIntersection::classifyHit
             // doing the surf2 with surf1 intersection but these
             // are merged later on)
 
-            label edge2I = getEdge(surf2, surf2Facei, nearLabel);
+            // edge hits all faces on surf2 connected to the edge
+            //
+            // The edge-edge intersection is symmetric, store only once.
+            // - When intersecting two surfaces, note which edges are cut each
+            //   time, but only create an edge from the first pass.
+            // - For self-intersection, it is slightly trickier if we don't
+            //   want too many duplicate points.
+
+            const label edge2I = getEdge(surf2, surf2Facei, nearLabel);
             const edge& e2 = surf2.edges()[edge2I];
+            label cutPointId = -1;
+
+            // Storage treatment
+            // =0: nothing/ignore
+            // >0: store point/edge-cut. Attempt to create new edge.
+            // <0: store point/edge-cut only
+            int handling = 0;
+            switch (cutFrom)
+            {
+                case surfaceIntersection::FIRST:
+                    handling = 1;
+                    break;
+                case surfaceIntersection::SECOND:
+                    handling = -1;
+                    break;
+                case surfaceIntersection::SELF:
+                    // The edge-edge intersection is hashed as an 'edge' to
+                    // exploit the commutative lookup.
+                    // Ie, only do the cut once
+                    const edge intersect(edgeI, edge2I);
+
+                    if (edgeEdgeIntersection_.insert(intersect))
+                    {
+                        handling = 1;
+                        forAll(e, edgepti)
+                        {
+                            const label endId = e[edgepti];
+                            const point& nearPt = surf1Pts[endId];
+
+                            if
+                            (
+                                mag(pHit.hitPoint() - nearPt)
+                              < surf1PointTol[endId]
+                            )
+                            {
+                                cutPointId = allCutPoints.size();
+
+                                if (avoidDuplicates_)
+                                {
+                                    if (edgeEndAsCut_.insert(endId, cutPointId))
+                                    {
+                                        // First time with this end-point
+                                        allCutPoints.append(nearPt);
+                                    }
+                                    else
+                                    {
+                                        // Already seen this end point
+                                        cutPointId = edgeEndAsCut_[endId];
+                                        handling = 2;  // cached
+                                    }
+                                }
+                                else
+                                {
+                                    allCutPoints.append(nearPt);
+                                }
+
+                                break;
+                            }
+                        }
+                    }
+                    break;
+            }
 
             if (debug & 2)
             {
-                Pout<< pHit.hitPoint() << " is surf1:"
-                    << " somewhere on edge " << e
-                    << " surf2: edge " << e2
-                    << " coords:" << surf2Pts[e2.start()]
-                    << surf2Pts[e2.end()] << endl;
+                Pout<< "hit-type[4] " << pHit.hitPoint() << " is surf1:"
+                    << " from edge[" << edgeI << "] " << e
+                    << "==" << e.line(surf1Pts)
+                    << " surf2: edge[" << edge2I << "] " << e2
+                    << " coords:" << e2.line(surf2Pts)
+                    << " - "
+                    << (
+                         handling < 0
+                       ? "cut-point" : handling
+                       ? "stored" : "suppressed"
+                       )
+                    << endl;
             }
 
-            allCutPoints.append(pHit.hitPoint());
-            surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
-
-            // edge hits all faces on surf2 connected to the edge
-
-            if (cutFrom == surfaceIntersection::FIRST)
+            if (handling)
             {
-                // edge-edge intersection is symmetric, store only
-                // once.
-                // edge hits all faces on surf2 connected to the
-                // edge
+                if (cutPointId == -1)
+                {
+                    cutPointId = allCutPoints.size();
+                    allCutPoints.append(pHit.hitPoint());
+                }
+                surfEdgeCuts[edgeI].append(cutPointId);
+            }
 
+            if (handling)
+            {
                 const labelList& facesB = surf2.edgeFaces()[edge2I];
-
                 forAll(facesB, faceBI)
                 {
                     storeIntersection
@@ -435,8 +557,9 @@ void Foam::surfaceIntersection::classifyHit
                         cutFrom,
                         facesA,
                         facesB[faceBI],
-                        allCutEdges,
-                        allCutPoints
+                        allCutPoints,
+                        cutPointId,
+                        allCutEdges
                     );
                 }
             }
@@ -448,77 +571,78 @@ void Foam::surfaceIntersection::classifyHit
         {
             // 5. Point hits face. Do what? Introduce
             // point & triangulation in face?
-            if (debug & 2)
-            {
-                Pout<< pHit.hitPoint() << " is surf1:"
-                    << " end point of edge " << e
-                    << " surf2: face " << surf2Facei
-                    << endl;
-            }
 
-            //
             // Look exactly at what side (of surf2) edge is. Leave out ones on
             // inside of surf2 (i.e. on opposite side of normal)
-            //
 
             // Vertex on/near surf2; vertex away from surf2
-            const label nearVert  = (edgeEnd == 0 ?  e.start() : e.end());
-            const label otherVert = (edgeEnd == 0 ?  e.end() : e.start());
+            // otherVert on outside of surf2
+            const label nearVert  = (edgeEnd == 0 ? e.start() : e.end());
+            const label otherVert = (edgeEnd == 0 ? e.end() : e.start());
 
-            const point& nearPt  = surf1.localPoints()[nearVert];
-            const point& otherPt = surf1.localPoints()[otherVert];
-
-            if (debug)
-            {
-                Pout
-                    << pHit.hitPoint() << " is surf1:"
-                    << " end point of edge " << e << " coord:"
-                    << nearPt
-                    << " surf2: face " << surf2Facei << endl;
-            }
+            const point& nearPt  = surf1Pts[nearVert];
+            const point& otherPt = surf1Pts[otherVert];
 
             const vector eVec = otherPt - nearPt;
 
             if ((surf2.faceNormals()[surf2Facei] & eVec) > 0)
             {
-                // otherVert on outside of surf2
+                // map to nearVert
+                // Reclassify as normal edge-face pierce (see below)
+                bool cached = false;
 
-                // Shift hitPoint a bit along edge.
-                //point hitPt = nearPt + 0.1*eVec;
-                point hitPt = nearPt;
+                label cutPointId = allCutPoints.size();
+                if (avoidDuplicates_)
+                {
+                    if (edgeEndAsCut_.insert(nearVert, cutPointId))
+                    {
+                        // First time with this end-point
+                        allCutPoints.append(nearPt);
+                    }
+                    else
+                    {
+                        // Already seen this end point
+                        cutPointId = edgeEndAsCut_[nearVert];
+                        cached = true;
+                    }
+                }
+                else
+                {
+                    allCutPoints.append(nearPt);
+                }
+
+                surfEdgeCuts[edgeI].append(cutPointId);
 
                 if (debug & 2)
                 {
-                    Pout<< "Shifted " << pHit.hitPoint()
-                        << " to " << hitPt
-                        << " along edge:" << e
-                        << " coords:" << surf1.localPoints()[e.start()]
-                        << surf1.localPoints()[e.end()] << endl;
+                    Pout<< "hit-type[5] " << pHit.hitPoint()
+                        << " shifted to " << nearPt
+                        << " from edge[" << edgeI << "] " << e
+                        << "==" << e.line(surf1Pts)
+                        << " hits surf2 face[" << surf2Facei << "]"
+                        << " - "
+                        << (cached ? "cached" : "stored") << endl;
                 }
 
-                // Reclassify as normal edge-face pierce (see below)
-
-                allCutPoints.append(hitPt);
-                surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
-
                 // edge hits single face only
                 storeIntersection
                 (
                     cutFrom,
                     facesA,
                     surf2Facei,
-                    allCutEdges,
-                    allCutPoints
+                    allCutPoints,
+                    cutPointId,
+                    allCutEdges
                 );
             }
             else
             {
                 if (debug & 2)
                 {
-                    Pout<< "Discarding " << pHit.hitPoint()
-                        << " since edge " << e << " on inside of surf2."
-                        << " surf2 normal:" << surf2.faceNormals()[surf2Facei]
-                        << endl;
+                    Pout<< "hit-type[5] " << pHit.hitPoint()
+                        << " from edge[" << edgeI << "] " << e
+                        << " hits inside of surf2 face[" << surf2Facei << "]"
+                        << " - discarded" << endl;
                 }
             }
         }
@@ -527,15 +651,16 @@ void Foam::surfaceIntersection::classifyHit
             // 6. Edge pierces face. 'Normal' situation.
             if (debug & 2)
             {
-                Pout<< pHit.hitPoint() << " is surf1:"
-                    << " somewhere on edge " << e
-                    << " surf2: face " << surf2Facei
-                    << endl;
+                Pout<< "hit-type[6] " << pHit.hitPoint()
+                    << " from edge[" << edgeI << "] " << e
+                    << "==" << e.line(surf1Pts)
+                    << " hits surf2 face[" << surf2Facei << "]"
+                    << " - stored" << endl;
             }
 
-            // edgeI intersects surf2. Store point.
+            const label cutPointId = allCutPoints.size();
             allCutPoints.append(pHit.hitPoint());
-            surfEdgeCuts[edgeI].append(allCutPoints.size()-1);
+            surfEdgeCuts[edgeI].append(cutPointId);
 
             // edge hits single face only
             storeIntersection
@@ -543,22 +668,19 @@ void Foam::surfaceIntersection::classifyHit
                 cutFrom,
                 facesA,
                 surf2Facei,
-                allCutEdges,
-                allCutPoints
+                allCutPoints,
+                cutPointId,
+                allCutEdges
             );
         }
     }
-    if (debug & 2)
-    {
-        Pout<< endl;
-    }
 }
 
 
 // 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 cutPoint
+// - facePairToVertex   : hash from face1I and face2I to (first) cutPoint
 // - facePairToEdge     : hash from face1I and face2I to cutEdge
 // - surfEdgeCuts       : gives for each edge the cutPoints
 //                        (in order from start to end)
@@ -569,12 +691,12 @@ void Foam::surfaceIntersection::doCutEdges
     const triSurfaceSearch& querySurf2,
     const enum originatingType cutFrom,
 
-    DynamicList<edge>& allCutEdges,
     DynamicList<point>& allCutPoints,
+    DynamicList<edge>& allCutEdges,
     List<DynamicList<label>>& surfEdgeCuts
 )
 {
-    const scalar oldTol = intersection::setPlanarTol(planarTol_);
+    const scalar oldTol = intersection::setPlanarTol(tolerance_);
 
     const pointField& surf1Pts = surf1.localPoints();
 
@@ -583,9 +705,7 @@ void Foam::surfaceIntersection::doCutEdges
 
     forAll(surf1PointTol, pointi)
     {
-        surf1PointTol[pointi] =
-            intersection::planarTol()
-          * minEdgeLen(surf1, pointi);
+        surf1PointTol[pointi] = tolerance_ * minEdgeLen(surf1, pointi);
     }
 
     const indexedOctree<treeDataPrimitivePatch<triSurface>>& searchTree
@@ -606,11 +726,11 @@ void Foam::surfaceIntersection::doCutEdges
             const edge& e = surf1.edges()[edgeI];
             const vector edgeVec = e.vec(surf1Pts);
 
-            // Extend start/end by tolerance - ensures cleaner cutting
+            // Extend start/end by 1/2 tolerance - ensures cleaner cutting
             const point ptStart =
-                surf1Pts[e.start()] - surf1PointTol[e.start()]*edgeVec;
+                surf1Pts[e.start()] - 0.5*surf1PointTol[e.start()]*edgeVec;
             const point ptEnd =
-                surf1Pts[e.end()]   + surf1PointTol[e.end()]*edgeVec;
+                surf1Pts[e.end()]   + 0.5*surf1PointTol[e.end()]*edgeVec;
 
             // Never intersect with faces attached to the edge itself
             maskFaces = surf1.edgeFaces()[edgeI];
@@ -639,8 +759,8 @@ void Foam::surfaceIntersection::doCutEdges
                     edgeI,
                     pHit,
 
-                    allCutEdges,
                     allCutPoints,
+                    allCutEdges,
                     surfEdgeCuts
                 );
             }
@@ -658,8 +778,11 @@ void Foam::surfaceIntersection::doCutEdges
             const point tolVec = intersection::planarTol()*(edgeVec);
             const scalar tolDim = mag(tolVec);
 
-            point ptStart = surf1Pts[e.start()];
-            const point ptEnd = surf1Pts[e.end()];
+            // Extend start/end by 1/2 tolerance - ensures cleaner cutting
+            point ptStart =
+                surf1Pts[e.start()] - 0.5*surf1PointTol[e.start()]*edgeVec;
+            const point ptEnd =
+                surf1Pts[e.end()]   + 0.5*surf1PointTol[e.end()]*edgeVec;
 
             bool doTrack = false;
             do
@@ -680,20 +803,21 @@ void Foam::surfaceIntersection::doCutEdges
                     edgeI,
                     pHit,
 
-                    allCutEdges,
                     allCutPoints,
+                    allCutEdges,
                     surfEdgeCuts
                 );
 
-                if (tolDim > 0.0)
+                if (tolerance_ > 0)
                 {
                     if (mag(pHit.hitPoint() - ptEnd) < tolDim)
                     {
+                        // Near the end => done
                         doTrack = false;
                     }
                     else
                     {
-                        // continue tracking a bit further on
+                        // Continue tracking a bit further on
                         ptStart = pHit.hitPoint() + tolVec;
                         doTrack = true;
                     }
@@ -702,6 +826,14 @@ void Foam::surfaceIntersection::doCutEdges
             while (doTrack);  // execute at least once
         }
     }
+    if (debug & 2)
+    {
+        Pout<< endl;
+    }
+
+    // These temporaries are now unneeded:
+    edgeEdgeIntersection_.clear();
+    edgeEndAsCut_.clear();
 
     intersection::setPlanarTol(oldTol);
 }
@@ -711,7 +843,10 @@ void Foam::surfaceIntersection::doCutEdges
 
 Foam::surfaceIntersection::surfaceIntersection()
 :
-    planarTol_(defaultTolerance),
+    tolerance_(1e-3),
+    allowEdgeHits_(true),
+    avoidDuplicates_(true),
+    warnDegenerate_(0),
     cutPoints_(0),
     cutEdges_(0),
     facePairToVertex_(0),
@@ -725,10 +860,13 @@ Foam::surfaceIntersection::surfaceIntersection
 (
     const triSurfaceSearch& query1,
     const triSurfaceSearch& query2,
-    const scalar planarTol
+    const dictionary& dict
 )
 :
-    planarTol_(planarTol),
+    tolerance_(1e-3),
+    allowEdgeHits_(true),
+    avoidDuplicates_(true),
+    warnDegenerate_(0),
     cutPoints_(0),
     cutEdges_(0),
     facePairToVertex_(2*max(query1.surface().size(), query2.surface().size())),
@@ -736,6 +874,8 @@ Foam::surfaceIntersection::surfaceIntersection
     surf1EdgeCuts_(0),
     surf2EdgeCuts_(0)
 {
+    setOptions(dict);
+
     const triSurface& surf1 = query1.surface();
     const triSurface& surf2 = query2.surface();
 
@@ -748,7 +888,7 @@ Foam::surfaceIntersection::surfaceIntersection
     }
 
 
-    DynamicList<edge> allCutEdges(surf1.nEdges()/20);
+    DynamicList<edge>  allCutEdges(surf1.nEdges()/20);
     DynamicList<point> allCutPoints(surf1.nPoints()/20);
 
 
@@ -761,8 +901,8 @@ Foam::surfaceIntersection::surfaceIntersection
         surf1,
         query2,
         surfaceIntersection::FIRST,
-        allCutEdges,
         allCutPoints,
+        allCutEdges,
         edgeCuts1
     );
     // Transfer to straight labelListList
@@ -786,8 +926,8 @@ Foam::surfaceIntersection::surfaceIntersection
         surf2,
         query1,
         surfaceIntersection::SECOND,
-        allCutEdges,
         allCutPoints,
+        allCutEdges,
         edgeCuts2
     );
 
@@ -818,16 +958,25 @@ Foam::surfaceIntersection::surfaceIntersection
         OFstream edge2Stream("surf2EdgeCuts.obj");
         writeIntersectedEdges(surf2, surf2EdgeCuts_, edge2Stream);
     }
+
+    // Temporaries
+    facePairToVertex_.clear();
+
+    // // Cleanup any duplicate cuts?
+    // mergeEdges();
 }
 
 
 Foam::surfaceIntersection::surfaceIntersection
 (
     const triSurfaceSearch& query1,
-    const scalar planarTol
+    const dictionary& dict
 )
 :
-    planarTol_(planarTol),
+    tolerance_(1e-3),
+    allowEdgeHits_(true),
+    avoidDuplicates_(true),
+    warnDegenerate_(0),
     cutPoints_(0),
     cutEdges_(0),
     facePairToVertex_(2*query1.surface().size()),
@@ -835,6 +984,8 @@ Foam::surfaceIntersection::surfaceIntersection
     surf1EdgeCuts_(0),
     surf2EdgeCuts_(0)
 {
+    setOptions(dict);
+
     const triSurface& surf1 = query1.surface();
 
     //
@@ -845,7 +996,7 @@ Foam::surfaceIntersection::surfaceIntersection
         Pout<< "Cutting surf1 edges" << endl;
     }
 
-    DynamicList<edge> allCutEdges;
+    DynamicList<edge>  allCutEdges;
     DynamicList<point> allCutPoints;
 
     // From edge to cut index on surface1
@@ -857,8 +1008,8 @@ Foam::surfaceIntersection::surfaceIntersection
         surf1,
         query1,
         surfaceIntersection::SELF,
-        allCutEdges,
         allCutPoints,
+        allCutEdges,
         edgeCuts1
     );
 
@@ -867,7 +1018,7 @@ Foam::surfaceIntersection::surfaceIntersection
     cutEdges_.transfer(allCutEdges);
     cutPoints_.transfer(allCutPoints);
 
-    // Shortcut.
+    // Short-circuit
     if (cutPoints_.empty() && cutEdges_.empty())
     {
         if (debug)
@@ -896,6 +1047,12 @@ Foam::surfaceIntersection::surfaceIntersection
         OFstream edge1Stream("surf1EdgeCuts.obj");
         writeIntersectedEdges(surf1, surf1EdgeCuts_, edge1Stream);
     }
+
+    // Temporaries
+    facePairToVertex_.clear();
+
+    // // Cleanup any duplicate cuts?
+    // mergeEdges();
 }
 
 
@@ -907,6 +1064,10 @@ Foam::surfaceIntersection::surfaceIntersection
     const edgeIntersections& intersections2
 )
 :
+    tolerance_(1e-3),
+    allowEdgeHits_(true),
+    avoidDuplicates_(true),
+    warnDegenerate_(0),
     cutPoints_(0),
     cutEdges_(0),
     facePairToVertex_(2*max(surf1.size(), surf2.size())),
@@ -916,7 +1077,7 @@ Foam::surfaceIntersection::surfaceIntersection
 {
 
     // All intersection Pout (so for both surfaces)
-    DynamicList<edge> allCutEdges((surf1.nEdges() + surf2.nEdges())/20);
+    DynamicList<edge>  allCutEdges((surf1.nEdges() + surf2.nEdges())/20);
     DynamicList<point> allCutPoints((surf1.nPoints() + surf2.nPoints())/20);
 
 
@@ -938,19 +1099,21 @@ Foam::surfaceIntersection::surfaceIntersection
 
             forAll(intersections, i)
             {
-                const pointIndexHit& pHit = intersections[i];
-
                 // edgeI intersects surf2. Store point.
+                const pointIndexHit& pHit = intersections[i];
+                const label cutPointId = allCutPoints.size();
+
                 allCutPoints.append(pHit.hitPoint());
-                edgeCuts1[edgeI].append(allCutPoints.size()-1);
+                edgeCuts1[edgeI].append(cutPointId);
 
                 storeIntersection
                 (
                     surfaceIntersection::FIRST,
                     surf1.edgeFaces()[edgeI],
-                    pHit.index(),               // surf2Facei
-                    allCutEdges,
-                    allCutPoints
+                    pHit.index(),
+                    allCutPoints,
+                    cutPointId,
+                    allCutEdges
                 );
             }
         }
@@ -960,7 +1123,6 @@ Foam::surfaceIntersection::surfaceIntersection
     }
 
 
-
     // Cut all edges of surf2 with surf1
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -979,19 +1141,21 @@ Foam::surfaceIntersection::surfaceIntersection
 
             forAll(intersections, i)
             {
-                const pointIndexHit& pHit = intersections[i];
-
                 // edgeI intersects surf1. Store point.
+                const pointIndexHit& pHit = intersections[i];
+                const label cutPointId = allCutPoints.size();
+
                 allCutPoints.append(pHit.hitPoint());
-                edgeCuts2[edgeI].append(allCutPoints.size()-1);
+                edgeCuts2[edgeI].append(cutPointId);
 
                 storeIntersection
                 (
                     surfaceIntersection::SECOND,
                     surf2.edgeFaces()[edgeI],
-                    pHit.index(),               // surf2Facei
-                    allCutEdges,
-                    allCutPoints
+                    pHit.index(),
+                    allCutPoints,
+                    cutPointId,
+                    allCutEdges
                 );
             }
         }
@@ -1037,8 +1201,7 @@ Foam::surfaceIntersection::surfaceIntersection
 
         forAllConstIter(labelPairLookup, facePairToEdge_, iter)
         {
-            label edgeI = iter();
-
+            const label edgeI = iter();
             const edge& e = cutEdges_[edgeI];
 
             usedPoints.insert(e[0]);
@@ -1047,7 +1210,7 @@ Foam::surfaceIntersection::surfaceIntersection
 
         forAllConstIter(labelPairLookup, facePairToVertex_, iter)
         {
-            label pointi = iter();
+            const label pointi = iter();
 
             if (!usedPoints.found(pointi))
             {
@@ -1058,6 +1221,12 @@ Foam::surfaceIntersection::surfaceIntersection
             }
         }
     }
+
+    // Temporaries
+    facePairToVertex_.clear();
+
+    // // Cleanup any duplicate cuts?
+    // mergeEdges();
 }
 
 
@@ -1075,12 +1244,6 @@ 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_;
@@ -1115,4 +1278,85 @@ const Foam::labelListList& Foam::surfaceIntersection::surf2EdgeCuts() const
 }
 
 
+void Foam::surfaceIntersection::mergePoints(const scalar mergeDist)
+{
+    pointField newPoints;
+    labelList pointMap;
+
+    const bool hasMerged = Foam::mergePoints
+    (
+        cutPoints_,
+        mergeDist,
+        false,
+        pointMap,
+        newPoints,
+        vector::zero
+    );
+
+    if (hasMerged)
+    {
+        cutPoints_.transfer(newPoints);
+
+        forAll(cutEdges_, edgei)
+        {
+            edge& e = cutEdges_[edgei];
+
+            e[0] = pointMap[e[0]];
+            e[1] = pointMap[e[1]];
+        }
+
+        forAll(surf1EdgeCuts_, edgei)
+        {
+            inplaceRenumber(pointMap, surf1EdgeCuts_[edgei]);
+            inplaceUniqueSort(surf1EdgeCuts_[edgei]);
+        }
+        forAll(surf2EdgeCuts_, edgei)
+        {
+            inplaceRenumber(pointMap, surf2EdgeCuts_[edgei]);
+            inplaceUniqueSort(surf2EdgeCuts_[edgei]);
+        }
+    }
+
+    this->mergeEdges();
+}
+
+
+void Foam::surfaceIntersection::mergeEdges()
+{
+    HashSet<edge, Hash<edge>> uniqEdges(2*cutEdges_.size());
+
+    label nUniqEdges = 0;
+    labelList edgeNumbering(cutEdges_.size(), -1);
+
+    forAll(cutEdges_, edgeI)
+    {
+        const edge& e = cutEdges_[edgeI];
+
+        // Remove degenerate and repeated edges
+        // - reordering (e[0] < e[1]) is not really necessary
+        if (e[0] != e[1] && uniqEdges.insert(e))
+        {
+            edgeNumbering[edgeI] = nUniqEdges;
+            if (nUniqEdges != edgeI)
+            {
+                cutEdges_[nUniqEdges] = e;
+            }
+            cutEdges_[nUniqEdges].sort();
+            ++nUniqEdges;
+        }
+    }
+
+    // if (nUniqEdges < cutEdges_.size())
+    // {
+    //     // Additional safety, in case the edge was replaced?
+    //     forAllIter(labelPairLookup, facePairToEdge_, iter)
+    //     {
+    //         iter() = edgeNumbering[iter()];
+    //     }
+    // }
+
+    cutEdges_.setSize(nUniqEdges);  // truncate
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
index ca8ba6ccc0..5c60958675 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.H
@@ -41,11 +41,19 @@ Description
     hit of both faces and an edge is created between the retrieved vertex and
     the new one.
 
-    Note: when doing intersecting itself uses intersection::planarTol() as a
-    fraction of
+    Note: when doing intersecting itself uses 'tolerance' as a fraction of
     current edge length to determine if intersection is a point-touching one
     instead of an edge-piercing action.
 
+    Some constructors allow a dictionary of intersection controls:
+    \table
+        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
+        warnDegenerate | Number of warnings about degenerate edges | label | 0
+    \endtable
+
 SourceFiles
     surfaceIntersection.C
     surfaceIntersectionFuncs.C
@@ -75,7 +83,7 @@ class triSurface;
 class edgeIntersections;
 
 /*---------------------------------------------------------------------------*\
-                    Class surfaceIntersection Declaration
+                     Class surfaceIntersection Declaration
 \*---------------------------------------------------------------------------*/
 
 class surfaceIntersection
@@ -91,13 +99,22 @@ class surfaceIntersection
         };
 
         //- Tolerance for intersections
-        scalar planarTol_;
+        scalar tolerance_;
+
+        //- Allow edge-ends to cut another edge.
+        bool allowEdgeHits_;
+
+        //- Avoid creating duplicate cuts near edge ends
+        bool avoidDuplicates_;
+
+        //- Maximum number of warnings about degenerate edges
+        label warnDegenerate_;
 
         //- Newly introduced points.
         pointField cutPoints_;
 
-        //- Newly introduced edges (are on both surfaces). Reference into
-        //  cutPoints.
+        //- Newly introduced edges (are on both surfaces).
+        //  Reference into cutPoints.
         edgeList cutEdges_;
 
         //- From face on surf1 and face on surf2 to intersection point
@@ -116,18 +133,27 @@ class surfaceIntersection
         //  If multiple cuts:sorted from edge.start to edge.end
         labelListList surf2EdgeCuts_;
 
+        //- Temporary storage to manage edge-edge self-intersections.
+        HashSet<edge, Hash<edge>> edgeEdgeIntersection_;
+
+        //- Temporary storage to manage cuts/intersections from the edge ends
+        Map<label> edgeEndAsCut_;
+
 
     // Private Member Functions
 
-        //- Write point in obj format.
-        static void writeOBJ(const point& pt, Ostream& os);
+        //- 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>&,
-            const List<edge>&,
-            Ostream&
+            const List<point>& pts,
+            const List<edge>& edges,
+            Ostream& os
         );
 
         //- Transfer contents of List<DynamicList<..>> to List<List<..>>
@@ -149,9 +175,6 @@ class surfaceIntersection
         //  to new (-1 if element removed)
         static void removeDuplicates(const labelList& map, labelList& labels);
 
-        //- Apply map to elements of a labelList
-        static void inlineRemap(const labelList& map, labelList& elems);
-
         // Remove all duplicate and degenerate elements. Return unique elements
         // and map from old to new.
         static edgeList filterEdges(const edgeList&, labelList& map);
@@ -159,30 +182,6 @@ class surfaceIntersection
         //- Remove all duplicate elements.
         static labelList filterLabels(const labelList& elems, labelList& map);
 
-        //- Do some checks if edge and face (resulting from hit)
-        //  should not be considered. Returns true if can be discarded.
-        static bool excludeEdgeHit
-        (
-            const triSurface& surf,
-            const label edgeI,
-            const label facei,
-            const scalar tol
-        );
-
-        ////- Given edge (eStart - eEnd) and normal direction construct plane
-        ////  and intersect all edges of hitFace with it.
-        ////  Return the edge and coordinate of hit.
-        //static pointIndexHit faceEdgeIntersection
-        //(
-        //    const triSurface&,
-        //    const label hitFacei,
-        //
-        //    const vector& n,
-        //    const point& eStart,
-        //    const point& eEnd
-        //);
-
-
         //- Debugging: Dump intersected edges to stream
         void writeIntersectedEdges
         (
@@ -199,7 +198,7 @@ class surfaceIntersection
             const scalar endTol,
             const point& p,
             const edge& e,
-            const pointField& points
+            const UList<point>& points
         );
 
         //- Update reference between faceA and faceB. Updates facePairToVertex_
@@ -209,8 +208,9 @@ class surfaceIntersection
             const enum originatingType cutFrom,
             const labelList& facesA,
             const label faceB,
-            DynamicList<edge>&,
-            DynamicList<point>&
+            const UList<point>& allCutPoints,
+            const label cutPointId,
+            DynamicList<edge>& allCutEdges
         );
 
         //- Investigate pHit to whether is case of point hits point,
@@ -224,8 +224,8 @@ class surfaceIntersection
             const label edgeI,
             const pointIndexHit& pHit,
 
-            DynamicList<edge>& allCutEdges,
             DynamicList<point>& allCutPoints,
+            DynamicList<edge>& allCutEdges,
             List<DynamicList<label>>& surfEdgeCuts
         );
 
@@ -236,8 +236,8 @@ class surfaceIntersection
             const triSurfaceSearch& querySurf2,
             const enum originatingType cutFrom,
 
-            DynamicList<edge>& allCutEdges,
             DynamicList<point>& allCutPoints,
+            DynamicList<edge>& allCutEdges,
             List<DynamicList<label>>& surfEdgeCuts
         );
 
@@ -246,9 +246,6 @@ public:
 
     // Public Data, Declarations
 
-    //- The default planarTol for intersections.
-    static const scalar defaultTolerance;
-
     ClassName("surfaceIntersection");
 
 
@@ -263,7 +260,7 @@ public:
         (
             const triSurfaceSearch& querySurf1,
             const triSurfaceSearch& querySurf2,
-            const scalar planarTol = surfaceIntersection::defaultTolerance
+            const dictionary& dict = dictionary::null
         );
 
         //- Construct from self-intersections.
@@ -271,7 +268,7 @@ public:
         surfaceIntersection
         (
             const triSurfaceSearch& querySurf1,
-            const scalar planarTol = surfaceIntersection::defaultTolerance
+            const dictionary& dict = dictionary::null
         );
 
         //- Construct from precalculated intersection information.
@@ -288,21 +285,33 @@ public:
 
     // Member Functions
 
+        //- The list of cut points
         const pointField& cutPoints() const;
 
+        //- The list of created edges
         const edgeList& cutEdges() const;
 
-        const labelPairLookup& facePairToVertex() const;
-
+        //- Lookup of pairs of faces to created edges
         const labelPairLookup& facePairToEdge() const;
 
         //- Access either surf1EdgeCuts (isFirstSurface = true) or
         //  surf2EdgeCuts
         const labelListList& edgeCuts(const bool isFirstSurf) const;
 
+        //- List of cut points on edges of surface1
         const labelListList& surf1EdgeCuts() const;
 
+        //- List of cut points on edges of surface2
         const labelListList& surf2EdgeCuts() const;
+
+
+        //- Geometric merge points (points within mergeDist) prior to
+        //  automatically calling mergeEdges().
+        void mergePoints(const scalar mergeDist);
+
+        //- Merge duplicate edges
+        void mergeEdges();
+
 };
 
 
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
index cd3af9da0c..6e261a86cb 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersectionFuncs.C
@@ -24,19 +24,24 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "surfaceIntersection.H"
+#include "triSurface.H"
 #include "triSurfaceSearch.H"
 #include "labelPairHashes.H"
 #include "OFstream.H"
-#include "HashSet.H"
-#include "triSurface.H"
-#include "pointIndexHit.H"
-#include "meshTools.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::surfaceIntersection::writeOBJ(const point& pt, Ostream& os)
+void Foam::surfaceIntersection::writeOBJ
+(
+    const List<point>& pts,
+    Ostream& os
+)
 {
-    os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
+    forAll(pts, i)
+    {
+        const point& pt = pts[i];
+        os << "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << nl;
+    }
 }
 
 
@@ -47,15 +52,13 @@ void Foam::surfaceIntersection::writeOBJ
     Ostream& os
 )
 {
-    forAll(pts, i)
-    {
-        writeOBJ(pts[i], os);
-    }
+    writeOBJ(pts, os);
+
     forAll(edges, i)
     {
         const edge& e = edges[i];
 
-        os << "l " << e.start()+1 << ' ' << e.end()+1 << endl;
+        os << "l " << e.start()+1 << ' ' << e.end()+1 << nl;
     }
 }
 
@@ -162,20 +165,6 @@ void Foam::surfaceIntersection::removeDuplicates
 }
 
 
-// Remap.
-void Foam::surfaceIntersection::inlineRemap
-(
-    const labelList& map,
-    labelList& elems
-)
-{
-    forAll(elems, elemI)
-    {
-        elems[elemI] = map[elems[elemI]];
-    }
-}
-
-
 // Remove all duplicate and degenerate elements. Return unique elements and
 // map from old to new.
 Foam::edgeList Foam::surfaceIntersection::filterEdges
@@ -265,14 +254,8 @@ void Foam::surfaceIntersection::writeIntersectedEdges
     // Dump all points (surface followed by cutPoints)
     const pointField& pts = surf.localPoints();
 
-    forAll(pts, pointi)
-    {
-        writeOBJ(pts[pointi], os);
-    }
-    forAll(cutPoints(), cutPointi)
-    {
-        writeOBJ(cutPoints()[cutPointi], os);
-    }
+    writeOBJ(pts, os);
+    writeOBJ(cutPoints(), os);
 
     forAll(edgeCutVerts, edgeI)
     {
@@ -284,16 +267,16 @@ void Foam::surfaceIntersection::writeIntersectedEdges
 
             // Start of original edge to first extra point
             os  << "l " << e.start()+1 << ' '
-                << extraVerts[0] + surf.nPoints() + 1 << endl;
+                << extraVerts[0] + surf.nPoints() + 1 << nl;
 
             for (label i = 1; i < extraVerts.size(); i++)
             {
                 os  << "l " << extraVerts[i-1] + surf.nPoints() + 1  << ' '
-                    << extraVerts[i] + surf.nPoints() + 1 << endl;
+                    << extraVerts[i] + surf.nPoints() + 1 << nl;
             }
 
             os  << "l " << extraVerts.last() + surf.nPoints() + 1
-                << ' ' << e.end()+1 << endl;
+                << ' ' << e.end()+1 << nl;
         }
     }
 }
@@ -306,7 +289,7 @@ Foam::label Foam::surfaceIntersection::classify
     const scalar endTol,
     const point& p,
     const edge& e,
-    const pointField& points
+    const UList<point>& points
 )
 {
     if (mag(p - points[e.start()]) < startTol)