diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C
index 5f85fc260..dae322af5 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,10 @@ Foam::automatic::automatic
     featureProximityFile_(coeffsDict_.lookup("featureProximityFile")),
     readInternalCloseness_(Switch(coeffsDict_.lookup("internalCloseness"))),
     internalClosenessFile_(coeffsDict_.lookup("internalClosenessFile")),
+    internalClosenessCellSizeCoeff_
+    (
+        readScalar(coeffsDict_.lookup("internalClosenessCellSizeCoeff"))
+    ),
     curvatureCellSizeCoeff_
     (
         readScalar(coeffsDict_.lookup("curvatureCellSizeCoeff"))
@@ -203,7 +207,7 @@ Foam::tmp<Foam::triSurfacePointScalarField> Foam::automatic::load()
         Info<< indent
             << "Reading internal closeness: " << internalClosenessFile_ << endl;
 
-        triSurfaceScalarField internalCloseness
+        triSurfacePointScalarField internalClosenessPointField
         (
             IOobject
             (
@@ -219,17 +223,13 @@ Foam::tmp<Foam::triSurfacePointScalarField> Foam::automatic::load()
             true
         );
 
-        scalarField internalClosenessPointField
-        (
-            patchInterpolate.faceToPointInterpolate(internalCloseness)
-        );
-
         forAll(pointCellSize, pI)
         {
             pointCellSize[pI] =
                 min
                 (
-                    internalClosenessPointField[meshPointMap[pI]],
+                    (1.0/internalClosenessCellSizeCoeff_)
+                   *internalClosenessPointField[meshPointMap[pI]],
                     pointCellSize[pI]
                 );
         }
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.H b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.H
index 07f4e9442..577f0617a 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.H
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,6 +74,10 @@ private:
         const Switch readInternalCloseness_;
         const word internalClosenessFile_;
 
+        //- The internalCloseness values are multiplied by the inverse
+        //  of this value to get the cell size
+        const scalar internalClosenessCellSizeCoeff_;
+
         //- The curvature values are multiplied by the inverse of this value to
         //  get the cell size
         const scalar curvatureCellSizeCoeff_;
diff --git a/applications/utilities/surface/surfaceFeatureExtract/Make/files b/applications/utilities/surface/surfaceFeatureExtract/Make/files
index a57125cd6..f8ead594d 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/Make/files
+++ b/applications/utilities/surface/surfaceFeatureExtract/Make/files
@@ -1,3 +1,6 @@
+surfaceExtractCloseness.C
+surfaceExtractPointCloseness.C
+surfaceFeatureExtractUtilities.C
 surfaceFeatureExtract.C
 
 EXE = $(FOAM_APPBIN)/surfaceFeatureExtract
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractCloseness.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractCloseness.C
new file mode 100644
index 000000000..1c6fb0521
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractCloseness.C
@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeatureExtract.H"
+#include "Time.H"
+#include "triSurfaceMesh.H"
+#include "vtkSurfaceWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::extractCloseness
+(
+    const fileName &sFeatFileName,
+    const Time& runTime,
+    const triSurface &surf,
+    const bool writeVTK
+)
+{
+    // Searchable triSurface
+    const triSurfaceMesh searchSurf
+    (
+        IOobject
+        (
+            sFeatFileName + ".closeness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf
+    );
+
+
+    // Prepare start and end points for intersection tests
+
+    const vectorField& normals = searchSurf.faceNormals();
+
+    const scalar span = searchSurf.bounds().mag();
+
+    const pointField start(searchSurf.faceCentres() - span*normals);
+    const pointField end(searchSurf.faceCentres() + span*normals);
+    const pointField& faceCentres = searchSurf.faceCentres();
+
+    List<List<pointIndexHit>> allHitinfo;
+
+    // Find all intersections (in order)
+    searchSurf.findLineAll(start, end, allHitinfo);
+
+    scalarField internalCloseness(start.size(), GREAT);
+    scalarField externalCloseness(start.size(), GREAT);
+
+    forAll(allHitinfo, fi)
+    {
+        const List<pointIndexHit>& hitInfo = allHitinfo[fi];
+
+        processHit
+        (
+            internalCloseness[fi],
+            externalCloseness[fi],
+            fi,
+            surf,
+            start[fi],
+            faceCentres[fi],
+            end[fi],
+            normals[fi],
+            normals,
+            hitInfo
+        );
+    }
+
+    triSurfaceScalarField internalClosenessField
+    (
+        IOobject
+        (
+            sFeatFileName + ".internalCloseness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf,
+        dimLength,
+        internalCloseness
+    );
+
+    internalClosenessField.write();
+
+    triSurfaceScalarField externalClosenessField
+    (
+        IOobject
+        (
+            sFeatFileName + ".externalCloseness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf,
+        dimLength,
+        externalCloseness
+    );
+
+    externalClosenessField.write();
+
+    if (writeVTK)
+    {
+        const faceList faces(surf.faces());
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constantPath()/"triSurface",// outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "internalCloseness",                // fieldName
+            internalCloseness,
+            false,                              // isNodeValues
+            true                                // verbose
+        );
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constantPath()/"triSurface",// outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "externalCloseness",                // fieldName
+            externalCloseness,
+            false,                              // isNodeValues
+            true                                // verbose
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractPointCloseness.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractPointCloseness.C
new file mode 100644
index 000000000..a0b2bde71
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceExtractPointCloseness.C
@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeatureExtract.H"
+#include "Time.H"
+#include "triSurfaceMesh.H"
+#include "vtkSurfaceWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::processHit
+(
+    scalar& internalCloseness,
+    scalar& externalCloseness,
+    const label fi,
+    const triSurface& surf,
+    const point& start,
+    const point& p,
+    const point& end,
+    const vector& normal,
+    const vectorField& normals,
+    const List<pointIndexHit>& hitInfo
+)
+{
+    if (hitInfo.size() < 1)
+    {
+        drawHitProblem(fi, surf, start, p, end, hitInfo);
+    }
+    else if (hitInfo.size() == 1)
+    {
+        if (!hitInfo[0].hit())
+        {
+        }
+        else if (hitInfo[0].index() != fi)
+        {
+            drawHitProblem(fi, surf, start, p, end, hitInfo);
+        }
+    }
+    else
+    {
+        label ownHiti = -1;
+
+        forAll(hitInfo, hI)
+        {
+            // Find the hit on the triangle that launched the ray
+
+            if (hitInfo[hI].index() == fi)
+            {
+                ownHiti = hI;
+                break;
+            }
+        }
+
+        if (ownHiti < 0)
+        {
+            drawHitProblem(fi, surf, start, p, end, hitInfo);
+        }
+        else if (ownHiti == 0)
+        {
+            // There are no internal hits, the first hit is the
+            // closest external hit
+
+            if
+            (
+                (normal & normals[hitInfo[ownHiti + 1].index()])
+              < externalToleranceCosAngle
+            )
+            {
+                externalCloseness = min
+                (
+                    externalCloseness,
+                    mag(p - hitInfo[ownHiti + 1].hitPoint())
+                );
+            }
+        }
+        else if (ownHiti == hitInfo.size() - 1)
+        {
+            // There are no external hits, the last but one hit is
+            // the closest internal hit
+
+            if
+            (
+                (normal & normals[hitInfo[ownHiti - 1].index()])
+              < internalToleranceCosAngle
+            )
+            {
+                internalCloseness = min
+                (
+                    internalCloseness,
+                    mag(p - hitInfo[ownHiti - 1].hitPoint())
+                );
+            }
+        }
+        else
+        {
+            if
+            (
+                (normal & normals[hitInfo[ownHiti + 1].index()])
+              < externalToleranceCosAngle
+            )
+            {
+                externalCloseness = min
+                (
+                    externalCloseness,
+                    mag(p - hitInfo[ownHiti + 1].hitPoint())
+                );
+            }
+
+            if
+            (
+                (normal & normals[hitInfo[ownHiti - 1].index()])
+              < internalToleranceCosAngle
+            )
+            {
+                internalCloseness = min
+                (
+                    internalCloseness,
+                    mag(p - hitInfo[ownHiti - 1].hitPoint())
+                );
+            }
+        }
+    }
+}
+
+
+void Foam::extractPointCloseness
+(
+    const fileName &sFeatFileName,
+    const Time& runTime,
+    const triSurface &surf,
+    const bool writeVTK
+)
+{
+    // Searchable triSurface
+    const triSurfaceMesh searchSurf
+    (
+        IOobject
+        (
+            sFeatFileName + ".closeness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf
+    );
+
+
+    // Prepare start and end points for intersection tests
+
+    const pointField& points = searchSurf.points();
+    const labelList& meshPoints = searchSurf.meshPoints();
+    const pointField& faceCentres = searchSurf.faceCentres();
+    const vectorField& normals = searchSurf.faceNormals();
+    const labelListList& pointFaces = searchSurf.pointFaces();
+
+    const scalar span = searchSurf.bounds().mag();
+
+    label nPointFaces = 0;
+    forAll(pointFaces, pfi)
+    {
+        nPointFaces += pointFaces[pfi].size();
+    }
+
+    pointField facePoints(nPointFaces);
+    pointField start(nPointFaces);
+    pointField end(nPointFaces);
+
+    label i = 0;
+    forAll(points, pi)
+    {
+        forAll(pointFaces[pi], pfi)
+        {
+            const label fi = pointFaces[pi][pfi];
+
+            facePoints[i] = (0.9*points[meshPoints[pi]] + 0.1*faceCentres[fi]);
+            const vector& n = normals[fi];
+
+            start[i] = facePoints[i] - span*n;
+            end[i] = facePoints[i] + span*n;
+
+            i++;
+        }
+    }
+
+    List<List<pointIndexHit>> allHitinfo;
+
+    // Find all intersections (in order)
+    searchSurf.findLineAll(start, end, allHitinfo);
+
+    scalarField internalCloseness(points.size(), GREAT);
+    scalarField externalCloseness(points.size(), GREAT);
+
+    i = 0;
+    forAll(points, pi)
+    {
+        forAll(pointFaces[pi], pfi)
+        {
+            const label fi = pointFaces[pi][pfi];
+            const List<pointIndexHit>& hitInfo = allHitinfo[i];
+
+            processHit
+            (
+                internalCloseness[pi],
+                externalCloseness[pi],
+                fi,
+                surf,
+                start[i],
+                facePoints[i],
+                end[i],
+                normals[fi],
+                normals,
+                hitInfo
+            );
+
+            i++;
+        }
+    }
+
+    triSurfacePointScalarField internalClosenessPointField
+    (
+        IOobject
+        (
+            sFeatFileName + ".internalPointCloseness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf,
+        dimLength,
+        internalCloseness
+    );
+
+    internalClosenessPointField.write();
+
+    triSurfacePointScalarField externalClosenessPointField
+    (
+        IOobject
+        (
+            sFeatFileName + ".externalPointCloseness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf,
+        dimLength,
+        externalCloseness
+    );
+
+    externalClosenessPointField.write();
+
+    if (writeVTK)
+    {
+        const faceList faces(surf.faces());
+        const Map<label>& meshPointMap = surf.meshPointMap();
+
+        forAll(meshPointMap, pi)
+        {
+            internalCloseness[pi] =
+                internalClosenessPointField[meshPointMap[pi]];
+
+            externalCloseness[pi] =
+                externalClosenessPointField[meshPointMap[pi]];
+        }
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constantPath()/"triSurface",// outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "internalPointCloseness",           // fieldName
+            internalCloseness,
+            true,                               // isNodeValues
+            true                                // verbose
+        );
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constantPath()/"triSurface",// outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "externalPointCloseness",           // fieldName
+            externalCloseness,
+            true,                               // isNodeValues
+            true                                // verbose
+        );
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 7395cf970..f435f3c7c 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -21,916 +21,19 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Application
-    surfaceFeatureExtract
-
-Description
-    Extracts and writes surface features to file. All but the basic feature
-    extraction is WIP.
-
-    Curvature calculation is an implementation of the algorithm from:
-
-      "Estimating Curvatures and their Derivatives on Triangle Meshes"
-      by S. Rusinkiewicz
-
 \*---------------------------------------------------------------------------*/
 
+#include "surfaceFeatureExtract.H"
 #include "argList.H"
 #include "Time.H"
-#include "triSurface.H"
-#include "surfaceFeatures.H"
 #include "featureEdgeMesh.H"
-#include "extendedFeatureEdgeMesh.H"
-#include "treeBoundBox.H"
-#include "meshTools.H"
-#include "OFstream.H"
-#include "triSurfaceMesh.H"
 #include "vtkSurfaceWriter.H"
-#include "triSurfaceFields.H"
-#include "indexedOctree.H"
-#include "treeDataEdge.H"
-#include "unitConversion.H"
-#include "plane.H"
-#include "tensor2D.H"
-#include "symmTensor2D.H"
-#include "point.H"
-#include "triadField.H"
-#include "transform.H"
 #include "IOdictionary.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-scalar calcVertexNormalWeight
-(
-    const triFace& f,
-    const label pI,
-    const vector& fN,
-    const pointField& points
-)
-{
-    label index = findIndex(f, pI);
-
-    if (index == -1)
-    {
-        FatalErrorInFunction
-            << "Point not in face" << abort(FatalError);
-    }
-
-    const vector e1 = points[f[index]] - points[f[f.fcIndex(index)]];
-    const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
-
-    return mag(fN)/(magSqr(e1)*magSqr(e2) + VSMALL);
-}
-
-
-point randomPointInPlane(const plane& p)
-{
-    // Perturb base point
-    const point& refPt = p.refPoint();
-
-    // ax + by + cz + d = 0
-    const FixedList<scalar, 4>& planeCoeffs = p.planeCoeffs();
-
-    const scalar perturbX = refPt.x() + 1e-3;
-    const scalar perturbY = refPt.y() + 1e-3;
-    const scalar perturbZ = refPt.z() + 1e-3;
-
-    if (mag(planeCoeffs[2]) < SMALL)
-    {
-        if (mag(planeCoeffs[1]) < SMALL)
-        {
-            const scalar x =
-                -1.0
-                *(
-                     planeCoeffs[3]
-                   + planeCoeffs[1]*perturbY
-                   + planeCoeffs[2]*perturbZ
-                 )/planeCoeffs[0];
-
-            return point
-            (
-                x,
-                perturbY,
-                perturbZ
-            );
-        }
-
-        const scalar y =
-            -1.0
-            *(
-                 planeCoeffs[3]
-               + planeCoeffs[0]*perturbX
-               + planeCoeffs[2]*perturbZ
-             )/planeCoeffs[1];
-
-        return point
-        (
-            perturbX,
-            y,
-            perturbZ
-        );
-    }
-    else
-    {
-        const scalar z =
-            -1.0
-            *(
-                 planeCoeffs[3]
-               + planeCoeffs[0]*perturbX
-               + planeCoeffs[1]*perturbY
-             )/planeCoeffs[2];
-
-        return point
-        (
-            perturbX,
-            perturbY,
-            z
-        );
-    }
-}
-
-
-triadField calcVertexCoordSys
-(
-    const triSurface& surf,
-    const vectorField& pointNormals
-)
-{
-    const pointField& points = surf.points();
-    const Map<label>& meshPointMap = surf.meshPointMap();
-
-    triadField pointCoordSys(points.size());
-
-    forAll(points, pI)
-    {
-        const point& pt = points[pI];
-        const vector& normal = pointNormals[meshPointMap[pI]];
-
-        if (mag(normal) < SMALL)
-        {
-            pointCoordSys[meshPointMap[pI]] = triad::unset;
-            continue;
-        }
-
-        plane p(pt, normal);
-
-        // Pick random point in plane
-        vector dir1 = pt - randomPointInPlane(p);
-        dir1 /= mag(dir1);
-
-        vector dir2 = dir1 ^ normal;
-        dir2 /= mag(dir2);
-
-        pointCoordSys[meshPointMap[pI]] = triad(dir1, dir2, normal);
-    }
-
-    return pointCoordSys;
-}
-
-
-vectorField calcVertexNormals(const triSurface& surf)
-{
-    // Weighted average of normals of faces attached to the vertex
-    // Weight = fA / (mag(e0)^2 * mag(e1)^2);
-
-    Info<< "Calculating vertex normals" << endl;
-
-    vectorField pointNormals(surf.nPoints(), Zero);
-
-    const pointField& points = surf.points();
-    const labelListList& pointFaces = surf.pointFaces();
-    const labelList& meshPoints = surf.meshPoints();
-
-    forAll(pointFaces, pI)
-    {
-        const labelList& pFaces = pointFaces[pI];
-
-        forAll(pFaces, fI)
-        {
-            const label facei = pFaces[fI];
-            const triFace& f = surf[facei];
-
-            vector fN = f.normal(points);
-
-            scalar weight = calcVertexNormalWeight
-            (
-                f,
-                meshPoints[pI],
-                fN,
-                points
-            );
-
-            pointNormals[pI] += weight*fN;
-        }
-
-        pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL;
-    }
-
-    return pointNormals;
-}
-
-
-triSurfacePointScalarField calcCurvature
-(
-    const word& name,
-    const Time& runTime,
-    const triSurface& surf,
-    const vectorField& pointNormals,
-    const triadField& pointCoordSys
-)
-{
-    Info<< "Calculating face curvature" << endl;
-
-    const pointField& points = surf.points();
-    const labelList& meshPoints = surf.meshPoints();
-    const Map<label>& meshPointMap = surf.meshPointMap();
-
-    triSurfacePointScalarField curvaturePointField
-    (
-        IOobject
-        (
-            name + ".curvature",
-            runTime.constant(),
-            "triSurface",
-            runTime,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        surf,
-        dimLength,
-        scalarField(points.size(), 0.0)
-    );
-
-    List<symmTensor2D> pointFundamentalTensors
-    (
-        points.size(),
-        symmTensor2D::zero
-    );
-
-    scalarList accumulatedWeights(points.size(), 0.0);
-
-    forAll(surf, fI)
-    {
-        const triFace& f = surf[fI];
-        const edgeList fEdges = f.edges();
-
-        // Calculate the edge vectors and the normal differences
-        vectorField edgeVectors(f.size(), Zero);
-        vectorField normalDifferences(f.size(), Zero);
-
-        forAll(fEdges, feI)
-        {
-            const edge& e = fEdges[feI];
-
-            edgeVectors[feI] = e.vec(points);
-            normalDifferences[feI] =
-               pointNormals[meshPointMap[e[0]]]
-             - pointNormals[meshPointMap[e[1]]];
-        }
-
-        // Set up a local coordinate system for the face
-        const vector& e0 = edgeVectors[0];
-        const vector eN = f.normal(points);
-        const vector e1 = (e0 ^ eN);
-
-        if (magSqr(eN) < ROOTVSMALL)
-        {
-            continue;
-        }
-
-        triad faceCoordSys(e0, e1, eN);
-        faceCoordSys.normalize();
-
-        // Construct the matrix to solve
-        scalarSymmetricSquareMatrix T(3, 0);
-        scalarDiagonalMatrix Z(3, 0);
-
-        // Least Squares
-        for (label i = 0; i < 3; ++i)
-        {
-            scalar x = edgeVectors[i] & faceCoordSys[0];
-            scalar y = edgeVectors[i] & faceCoordSys[1];
-
-            T(0, 0) += sqr(x);
-            T(1, 0) += x*y;
-            T(1, 1) += sqr(x) + sqr(y);
-            T(2, 1) += x*y;
-            T(2, 2) += sqr(y);
-
-            scalar dndx = normalDifferences[i] & faceCoordSys[0];
-            scalar dndy = normalDifferences[i] & faceCoordSys[1];
-
-            Z[0] += dndx*x;
-            Z[1] += dndx*y + dndy*x;
-            Z[2] += dndy*y;
-        }
-
-        // Perform Cholesky decomposition and back substitution.
-        // Decomposed matrix is in T and solution is in Z.
-        LUsolve(T, Z);
-        symmTensor2D secondFundamentalTensor(Z[0], Z[1], Z[2]);
-
-        // Loop over the face points adding the contribution of the face
-        // curvature to the points.
-        forAll(f, fpI)
-        {
-            const label patchPointIndex = meshPointMap[f[fpI]];
-
-            const triad& ptCoordSys = pointCoordSys[patchPointIndex];
-
-            if (!ptCoordSys.set())
-            {
-                continue;
-            }
-
-            // Rotate faceCoordSys to ptCoordSys
-            tensor rotTensor = rotationTensor(ptCoordSys[2], faceCoordSys[2]);
-            triad rotatedFaceCoordSys = rotTensor & tensor(faceCoordSys);
-
-            // Project the face curvature onto the point plane
-
-            vector2D cmp1
-            (
-                ptCoordSys[0] & rotatedFaceCoordSys[0],
-                ptCoordSys[0] & rotatedFaceCoordSys[1]
-            );
-            vector2D cmp2
-            (
-                ptCoordSys[1] & rotatedFaceCoordSys[0],
-                ptCoordSys[1] & rotatedFaceCoordSys[1]
-            );
-
-            tensor2D projTensor
-            (
-                cmp1,
-                cmp2
-            );
-
-            symmTensor2D projectedFundamentalTensor
-            (
-                projTensor.x() & (secondFundamentalTensor & projTensor.x()),
-                projTensor.x() & (secondFundamentalTensor & projTensor.y()),
-                projTensor.y() & (secondFundamentalTensor & projTensor.y())
-            );
-
-            // Calculate weight
-            // TODO: Voronoi area weighting
-            scalar weight = calcVertexNormalWeight
-            (
-                f,
-                meshPoints[patchPointIndex],
-                f.normal(points),
-                points
-            );
-
-            // Sum contribution of face to this point
-            pointFundamentalTensors[patchPointIndex] +=
-                weight*projectedFundamentalTensor;
-
-            accumulatedWeights[patchPointIndex] += weight;
-        }
-
-        if (false)
-        {
-            Info<< "Points = " << points[f[0]] << " "
-                << points[f[1]] << " "
-                << points[f[2]] << endl;
-            Info<< "edgeVecs = " << edgeVectors[0] << " "
-                << edgeVectors[1] << " "
-                << edgeVectors[2] << endl;
-            Info<< "normDiff = " << normalDifferences[0] << " "
-                << normalDifferences[1] << " "
-                << normalDifferences[2] << endl;
-            Info<< "faceCoordSys = " << faceCoordSys << endl;
-            Info<< "T = " << T << endl;
-            Info<< "Z = " << Z << endl;
-        }
-    }
-
-    forAll(curvaturePointField, pI)
-    {
-        pointFundamentalTensors[pI] /= (accumulatedWeights[pI] + SMALL);
-
-        vector2D principalCurvatures = eigenValues(pointFundamentalTensors[pI]);
-
-        //scalar curvature =
-        //    (principalCurvatures[0] + principalCurvatures[1])/2;
-        scalar curvature = max
-        (
-            mag(principalCurvatures[0]),
-            mag(principalCurvatures[1])
-        );
-        //scalar curvature = principalCurvatures[0]*principalCurvatures[1];
-
-        curvaturePointField[meshPoints[pI]] = curvature;
-    }
-
-    return curvaturePointField;
-}
-
-
-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,
-    const scalar defaultCellSize
-)
-{
-    scalar minDist = defaultCellSize;
-
-    for
-    (
-        label hI1 = 0;
-        hI1 < hitList.size() - 1;
-        ++hI1
-    )
-    {
-        const pointIndexHit& pHit1 = hitList[hI1];
-
-        if (pHit1.hit())
-        {
-            for
-            (
-                label hI2 = hI1 + 1;
-                hI2 < hitList.size();
-                ++hI2
-            )
-            {
-                const pointIndexHit& pHit2 = hitList[hI2];
-
-                if (pHit2.hit())
-                {
-                    scalar curDist = mag(pHit1.hitPoint() - pHit2.hitPoint());
-
-                    minDist = min(curDist, minDist);
-                }
-            }
-        }
-    }
-
-    return minDist;
-}
-
-
-scalar calcProximityOfFeatureEdges
-(
-    const extendedFeatureEdgeMesh& efem,
-    const List<pointIndexHit>& hitList,
-    const scalar defaultCellSize
-)
-{
-    scalar minDist = defaultCellSize;
-
-    for
-    (
-        label hI1 = 0;
-        hI1 < hitList.size() - 1;
-        ++hI1
-    )
-    {
-        const pointIndexHit& pHit1 = hitList[hI1];
-
-        if (pHit1.hit())
-        {
-            const edge& e1 = efem.edges()[pHit1.index()];
-
-            for
-            (
-                label hI2 = hI1 + 1;
-                hI2 < hitList.size();
-                ++hI2
-            )
-            {
-                const pointIndexHit& pHit2 = hitList[hI2];
-
-                if (pHit2.hit())
-                {
-                    const edge& e2 = efem.edges()[pHit2.index()];
-
-                    // Don't refine if the edges are connected to each other
-                    if (!edgesConnected(e1, e2))
-                    {
-                        scalar curDist =
-                            mag(pHit1.hitPoint() - pHit2.hitPoint());
-
-                        minDist = min(curDist, minDist);
-                    }
-                }
-            }
-        }
-    }
-
-    return minDist;
-}
-
-
-void dumpBox(const treeBoundBox& bb, const fileName& fName)
-{
-    OFstream str(fName);
-
-    Info<< "Dumping bounding box " << bb << " as lines to obj file "
-        << str.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;
-    }
-}
-
-
-// Deletes all edges inside/outside bounding box from set.
-void deleteBox
-(
-    const triSurface& surf,
-    const treeBoundBox& bb,
-    const bool removeInside,
-    List<surfaceFeatures::edgeStatus>& edgeStat
-)
-{
-    forAll(edgeStat, edgeI)
-    {
-        const point eMid = surf.edges()[edgeI].centre(surf.localPoints());
-
-        if (removeInside ? bb.contains(eMid) : !bb.contains(eMid))
-        {
-            edgeStat[edgeI] = surfaceFeatures::NONE;
-        }
-    }
-}
-
-
-bool onLine(const point& p, const linePointRef& line)
-{
-    const point& a = line.start();
-    const point& b = line.end();
-
-    if
-    (
-        ( p.x() < min(a.x(), b.x()) || p.x() > max(a.x(), b.x()) )
-     || ( p.y() < min(a.y(), b.y()) || p.y() > max(a.y(), b.y()) )
-     || ( p.z() < min(a.z(), b.z()) || p.z() > max(a.z(), b.z()) )
-    )
-    {
-        return false;
-    }
-
-    return true;
-}
-
-
-// Deletes all edges inside/outside bounding box from set.
-void deleteEdges
-(
-    const triSurface& surf,
-    const plane& cutPlane,
-    List<surfaceFeatures::edgeStatus>& edgeStat
-)
-{
-    const pointField& points = surf.points();
-    const labelList& meshPoints = surf.meshPoints();
-
-    forAll(edgeStat, edgeI)
-    {
-        const edge& e = surf.edges()[edgeI];
-        const point& p0 = points[meshPoints[e.start()]];
-        const point& p1 = points[meshPoints[e.end()]];
-        const linePointRef line(p0, p1);
-
-        // If edge does not intersect the plane, delete.
-        scalar intersect = cutPlane.lineIntersect(line);
-
-        point featPoint = intersect * (p1 - p0) + p0;
-
-        if (!onLine(featPoint, line))
-        {
-            edgeStat[edgeI] = surfaceFeatures::NONE;
-        }
-    }
-}
-
-
-void drawHitProblem
-(
-    label fI,
-    const triSurface& surf,
-    const pointField& start,
-    const pointField& faceCentres,
-    const pointField& end,
-    const List<pointIndexHit>& hitInfo
-)
-{
-    Info<< nl << "# findLineAll did not hit its own face."
-        << nl << "# fI " << fI
-        << nl << "# start " << start[fI]
-        << nl << "# f centre " << faceCentres[fI]
-        << nl << "# end " << end[fI]
-        << nl << "# hitInfo " << hitInfo
-        << endl;
-
-    meshTools::writeOBJ(Info, start[fI]);
-    meshTools::writeOBJ(Info, faceCentres[fI]);
-    meshTools::writeOBJ(Info, end[fI]);
-
-    Info<< "l 1 2 3" << endl;
-
-    meshTools::writeOBJ(Info, surf.points()[surf[fI][0]]);
-    meshTools::writeOBJ(Info, surf.points()[surf[fI][1]]);
-    meshTools::writeOBJ(Info, surf.points()[surf[fI][2]]);
-
-    Info<< "f 4 5 6" << endl;
-
-    forAll(hitInfo, hI)
-    {
-        label hFI = hitInfo[hI].index();
-
-        meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
-        meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
-        meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
-
-        Info<< "f "
-            << 3*hI + 7 << " "
-            << 3*hI + 8 << " "
-            << 3*hI + 9
-            << endl;
-    }
-}
-
-
-// Unmark non-manifold edges if individual triangles are not features
-void unmarkBaffles
-(
-    const triSurface& surf,
-    const scalar includedAngle,
-    List<surfaceFeatures::edgeStatus>& edgeStat
-)
-{
-    scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
-
-    const labelListList& edgeFaces = surf.edgeFaces();
-
-    forAll(edgeFaces, edgeI)
-    {
-        const labelList& eFaces = edgeFaces[edgeI];
-
-        if (eFaces.size() > 2)
-        {
-            label i0 = eFaces[0];
-            //const labelledTri& f0 = surf[i0];
-            const Foam::vector& n0 = surf.faceNormals()[i0];
-
-            //Pout<< "edge:" << edgeI << " n0:" << n0 << endl;
-
-            bool same = true;
-
-            for (label i = 1; i < eFaces.size(); i++)
-            {
-                //const labelledTri& f = surf[i];
-                const Foam::vector& n = surf.faceNormals()[eFaces[i]];
-
-                //Pout<< "    mag(n&n0): " << mag(n&n0) << endl;
-
-                if (mag(n&n0) < minCos)
-                {
-                    same = false;
-                    break;
-                }
-            }
-
-            if (same)
-            {
-                edgeStat[edgeI] = surfaceFeatures::NONE;
-            }
-        }
-    }
-}
-
-
-//- Divide into multiple normal bins
-//  - return REGION if != 2 normals
-//  - return REGION if 2 normals that make feature angle
-//  - otherwise return NONE and set normals,bins
-surfaceFeatures::edgeStatus checkFlatRegionEdge
-(
-    const triSurface& surf,
-    const scalar tol,
-    const scalar includedAngle,
-    const label edgeI
-)
-{
-    const edge& e = surf.edges()[edgeI];
-    const labelList& eFaces = surf.edgeFaces()[edgeI];
-
-    // Bin according to normal
-
-    DynamicList<Foam::vector> normals(2);
-    DynamicList<labelList> bins(2);
-
-    forAll(eFaces, eFacei)
-    {
-        const Foam::vector& n = surf.faceNormals()[eFaces[eFacei]];
-
-        // Find the normal in normals
-        label index = -1;
-        forAll(normals, normalI)
-        {
-            if (mag(n&normals[normalI]) > (1-tol))
-            {
-                index = normalI;
-                break;
-            }
-        }
-
-        if (index != -1)
-        {
-            bins[index].append(eFacei);
-        }
-        else if (normals.size() >= 2)
-        {
-            // Would be third normal. Mark as feature.
-            //Pout<< "** at edge:" << surf.localPoints()[e[0]]
-            //    << surf.localPoints()[e[1]]
-            //    << " have normals:" << normals
-            //    << " and " << n << endl;
-            return surfaceFeatures::REGION;
-        }
-        else
-        {
-            normals.append(n);
-            bins.append(labelList(1, eFacei));
-        }
-    }
-
-
-    // Check resulting number of bins
-    if (bins.size() == 1)
-    {
-        // Note: should check here whether they are two sets of faces
-        // that are planar or indeed 4 faces al coming together at an edge.
-        //Pout<< "** at edge:"
-        //    << surf.localPoints()[e[0]]
-        //    << surf.localPoints()[e[1]]
-        //    << " have single normal:" << normals[0]
-        //    << endl;
-        return surfaceFeatures::NONE;
-    }
-    else
-    {
-        // Two bins. Check if normals make an angle
-
-        //Pout<< "** at edge:"
-        //    << surf.localPoints()[e[0]]
-        //    << surf.localPoints()[e[1]] << nl
-        //    << "    normals:" << normals << nl
-        //    << "    bins   :" << bins << nl
-        //    << endl;
-
-        if (includedAngle >= 0)
-        {
-            scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
-
-            forAll(eFaces, i)
-            {
-                const Foam::vector& ni = surf.faceNormals()[eFaces[i]];
-                for (label j=i+1; j<eFaces.size(); j++)
-                {
-                    const Foam::vector& nj = surf.faceNormals()[eFaces[j]];
-                    if (mag(ni & nj) < minCos)
-                    {
-                        //Pout<< "have sharp feature between normal:" << ni
-                        //    << " and " << nj << endl;
-
-                        // Is feature. Keep as region or convert to
-                        // feature angle? For now keep as region.
-                        return surfaceFeatures::REGION;
-                    }
-                }
-            }
-        }
-
-
-        // So now we have two normals bins but need to make sure both
-        // bins have the same regions in it.
-
-         // 1. store + or - region number depending
-        //    on orientation of triangle in bins[0]
-        const labelList& bin0 = bins[0];
-        labelList regionAndNormal(bin0.size());
-        forAll(bin0, i)
-        {
-            const labelledTri& t = surf.localFaces()[eFaces[bin0[i]]];
-            int dir = t.edgeDirection(e);
-
-            if (dir > 0)
-            {
-                regionAndNormal[i] = t.region()+1;
-            }
-            else if (dir == 0)
-            {
-                FatalErrorInFunction
-                    << exit(FatalError);
-            }
-            else
-            {
-                regionAndNormal[i] = -(t.region()+1);
-            }
-        }
-
-        // 2. check against bin1
-        const labelList& bin1 = bins[1];
-        labelList regionAndNormal1(bin1.size());
-        forAll(bin1, i)
-        {
-            const labelledTri& t = surf.localFaces()[eFaces[bin1[i]]];
-            int dir = t.edgeDirection(e);
-
-            label myRegionAndNormal;
-            if (dir > 0)
-            {
-                myRegionAndNormal = t.region()+1;
-            }
-            else
-            {
-                myRegionAndNormal = -(t.region()+1);
-            }
-
-            regionAndNormal1[i] = myRegionAndNormal;
-
-            label index = findIndex(regionAndNormal, -myRegionAndNormal);
-            if (index == -1)
-            {
-                // Not found.
-                //Pout<< "cannot find region " << myRegionAndNormal
-                //    << " in regions " << regionAndNormal << endl;
-
-                return surfaceFeatures::REGION;
-            }
-        }
-
-        // Passed all checks, two normal bins with the same contents.
-        //Pout<< "regionAndNormal:" << regionAndNormal << endl;
-        //Pout<< "myRegionAndNormal:" << regionAndNormal1 << endl;
-
-        return surfaceFeatures::NONE;
-    }
-}
-
-
-void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
-{
-    os  << "    points : " << fem.points().size() << nl
-        << "    of which" << nl
-        << "        convex             : "
-        << fem.concaveStart() << nl
-        << "        concave            : "
-        << (fem.mixedStart()-fem.concaveStart()) << nl
-        << "        mixed              : "
-        << (fem.nonFeatureStart()-fem.mixedStart()) << nl
-        << "        non-feature        : "
-        << (fem.points().size()-fem.nonFeatureStart()) << nl
-        << "    edges  : " << fem.edges().size() << nl
-        << "    of which" << nl
-        << "        external edges     : "
-        << fem.internalStart() << nl
-        << "        internal edges     : "
-        << (fem.flatStart()- fem.internalStart()) << nl
-        << "        flat edges         : "
-        << (fem.openStart()- fem.flatStart()) << nl
-        << "        open edges         : "
-        << (fem.multipleStart()- fem.openStart()) << nl
-        << "        multiply connected : "
-        << (fem.edges().size()- fem.multipleStart()) << endl;
-}
-
-
-
 int main(int argc, char *argv[])
 {
     argList::addNote
@@ -997,13 +100,7 @@ int main(int argc, char *argv[])
         surf.writeStats(Info);
         Info<< endl;
 
-        faceList faces(surf.size());
-
-        forAll(surf, fI)
-        {
-            faces[fI] = surf[fI].triFaceFace();
-        }
-
+        const faceList faces(surf.faces());
 
         // Either construct features from surface & featureAngle or read set.
         // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1224,11 +321,6 @@ int main(int argc, char *argv[])
             << "        internal edges : " << newSet.nInternalEdges() << nl
             << endl;
 
-        //if (writeObj)
-        //{
-        //    newSet.writeObj("final");
-        //}
-
         boolList surfBaffleRegions(surf.patches().size(), false);
 
         wordList surfBaffleNames;
@@ -1318,349 +410,20 @@ int main(int argc, char *argv[])
 
         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
-    //     )
-    // );
 
+        // Find distance between close features
         if (closeness)
         {
             Info<< nl << "Extracting internal and external closeness of "
                 << "surface." << endl;
 
-
-            triSurfaceMesh searchSurf
-            (
-                IOobject
-                (
-                    sFeatFileName + ".closeness",
-                    runTime.constant(),
-                    "triSurface",
-                    runTime,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                surf
-            );
-
-
-            // Internal and external closeness
-
-            // Prepare start and end points for intersection tests
-
-            const vectorField& normals = searchSurf.faceNormals();
-
-            scalar span = searchSurf.bounds().mag();
-
-            scalar externalAngleTolerance = 10;
-            scalar externalToleranceCosAngle =
-                Foam::cos
-                (
-                    degToRad(180 - externalAngleTolerance)
-                );
-
-            scalar internalAngleTolerance = 45;
-            scalar internalToleranceCosAngle =
-                Foam::cos
-                (
-                    degToRad(180 - internalAngleTolerance)
-                );
-
             Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle
                 << nl
                 << "internalToleranceCosAngle: " << internalToleranceCosAngle
                 << endl;
 
-            // Info<< "span " << span << endl;
-
-            pointField start(searchSurf.faceCentres() - span*normals);
-            pointField end(searchSurf.faceCentres() + span*normals);
-            const pointField& faceCentres = searchSurf.faceCentres();
-
-            List<List<pointIndexHit>> allHitInfo;
-
-            // Find all intersections (in order)
-            searchSurf.findLineAll(start, end, allHitInfo);
-
-            scalarField internalCloseness(start.size(), GREAT);
-            scalarField externalCloseness(start.size(), GREAT);
-
-            forAll(allHitInfo, fI)
-            {
-                const List<pointIndexHit>& hitInfo = allHitInfo[fI];
-
-                if (hitInfo.size() < 1)
-                {
-                    drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
-
-                    // FatalErrorIn(args.executable())
-                    //     << "findLineAll did not hit its own face."
-                    //     << exit(FatalError);
-                }
-                else if (hitInfo.size() == 1)
-                {
-                    if (!hitInfo[0].hit())
-                    {
-                        // FatalErrorIn(args.executable())
-                        //     << "findLineAll did not hit any face."
-                        //     << exit(FatalError);
-                    }
-                    else if (hitInfo[0].index() != fI)
-                    {
-                        drawHitProblem
-                        (
-                            fI,
-                            surf,
-                            start,
-                            faceCentres,
-                            end,
-                            hitInfo
-                        );
-
-                        // FatalErrorIn(args.executable())
-                        //     << "findLineAll did not hit its own face."
-                        //     << exit(FatalError);
-                    }
-                }
-                else
-                {
-                    label ownHitI = -1;
-
-                    forAll(hitInfo, hI)
-                    {
-                        // Find the hit on the triangle that launched the ray
-
-                        if (hitInfo[hI].index() == fI)
-                        {
-                            ownHitI = hI;
-
-                            break;
-                        }
-                    }
-
-                    if (ownHitI < 0)
-                    {
-                        drawHitProblem
-                        (
-                            fI,
-                            surf,
-                            start,
-                            faceCentres,
-                            end,
-                            hitInfo
-                        );
-
-                        // FatalErrorIn(args.executable())
-                        //     << "findLineAll did not hit its own face."
-                        //     << exit(FatalError);
-                    }
-                    else if (ownHitI == 0)
-                    {
-                        // There are no internal hits, the first hit is the
-                        // closest external hit
-
-                        if
-                        (
-                            (
-                                normals[fI]
-                              & normals[hitInfo[ownHitI + 1].index()]
-                            )
-                          < externalToleranceCosAngle
-                        )
-                        {
-                            externalCloseness[fI] =
-                                mag
-                                (
-                                    faceCentres[fI]
-                                  - hitInfo[ownHitI + 1].hitPoint()
-                                );
-                        }
-                    }
-                    else if (ownHitI == hitInfo.size() - 1)
-                    {
-                        // There are no external hits, the last but one hit is
-                        // the closest internal hit
-
-                        if
-                        (
-                            (
-                                normals[fI]
-                              & normals[hitInfo[ownHitI - 1].index()]
-                            )
-                          < internalToleranceCosAngle
-                        )
-                        {
-                            internalCloseness[fI] =
-                                mag
-                                (
-                                    faceCentres[fI]
-                                  - hitInfo[ownHitI - 1].hitPoint()
-                                );
-                        }
-                    }
-                    else
-                    {
-                        if
-                        (
-                            (
-                                normals[fI]
-                              & normals[hitInfo[ownHitI + 1].index()]
-                            )
-                          < externalToleranceCosAngle
-                        )
-                        {
-                            externalCloseness[fI] =
-                                mag
-                                (
-                                    faceCentres[fI]
-                                  - hitInfo[ownHitI + 1].hitPoint()
-                                );
-                        }
-
-                        if
-                        (
-                            (
-                                normals[fI]
-                              & normals[hitInfo[ownHitI - 1].index()]
-                            )
-                          < internalToleranceCosAngle
-                        )
-                        {
-                            internalCloseness[fI] =
-                                mag
-                                (
-                                    faceCentres[fI]
-                                  - hitInfo[ownHitI - 1].hitPoint()
-                                );
-                        }
-                    }
-                }
-            }
-
-            triSurfaceScalarField internalClosenessField
-            (
-                IOobject
-                (
-                    sFeatFileName + ".internalCloseness",
-                    runTime.constant(),
-                    "triSurface",
-                    runTime,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                surf,
-                dimLength,
-                internalCloseness
-            );
-
-            internalClosenessField.write();
-
-            triSurfaceScalarField externalClosenessField
-            (
-                IOobject
-                (
-                    sFeatFileName + ".externalCloseness",
-                    runTime.constant(),
-                    "triSurface",
-                    runTime,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                surf,
-                dimLength,
-                externalCloseness
-            );
-
-            externalClosenessField.write();
-
-            if (writeVTK)
-            {
-                vtkSurfaceWriter().write
-                (
-                    runTime.constantPath()/"triSurface",// outputDir
-                    sFeatFileName,                      // surfaceName
-                    surf.points(),
-                    faces,
-                    "internalCloseness",                // fieldName
-                    internalCloseness,
-                    false,                              // isNodeValues
-                    true                                // verbose
-                );
-
-                vtkSurfaceWriter().write
-                (
-                    runTime.constantPath()/"triSurface",// outputDir
-                    sFeatFileName,                      // surfaceName
-                    surf.points(),
-                    faces,
-                    "externalCloseness",                // fieldName
-                    externalCloseness,
-                    false,                              // isNodeValues
-                    true                                // verbose
-                );
-            }
+            extractCloseness(sFeatFileName, runTime, surf, writeVTK);
+            extractPointCloseness(sFeatFileName, runTime, surf, writeVTK);
         }
 
 
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H
new file mode 100644
index 000000000..38701a9dc
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H
@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    surfaceFeatureExtract
+
+Description
+    Utility functions for surfaceFeatureExtract
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeatures.H"
+#include "extendedFeatureEdgeMesh.H"
+#include "triSurfaceFields.H"
+#include "plane.H"
+#include "triadField.H"
+
+namespace Foam
+{
+    extern const scalar internalAngleTolerance;
+    extern const scalar internalToleranceCosAngle;
+
+    extern const scalar externalAngleTolerance;
+    extern const scalar externalToleranceCosAngle;
+
+    scalar calcVertexNormalWeight
+    (
+        const triFace& f,
+        const label pI,
+        const vector& fN,
+        const pointField& points
+    );
+
+    point randomPointInPlane(const plane& p);
+
+    triadField calcVertexCoordSys
+    (
+        const triSurface& surf,
+        const vectorField& pointNormals
+    );
+
+    vectorField calcVertexNormals(const triSurface& surf);
+
+    triSurfacePointScalarField calcCurvature
+    (
+        const word& name,
+        const Time& runTime,
+        const triSurface& surf,
+        const vectorField& pointNormals,
+        const triadField& pointCoordSys
+    );
+
+    bool edgesConnected(const edge& e1, const edge& e2);
+
+    scalar calcProximityOfFeaturePoints
+    (
+        const List<pointIndexHit>& hitList,
+        const scalar defaultCellSize
+    );
+
+    scalar calcProximityOfFeatureEdges
+    (
+        const extendedFeatureEdgeMesh& efem,
+        const List<pointIndexHit>& hitList,
+        const scalar defaultCellSize
+    );
+
+    void dumpBox(const treeBoundBox& bb, const fileName& fName);
+
+    //- Deletes all edges inside/outside bounding box from set.
+    void deleteBox
+    (
+        const triSurface& surf,
+        const treeBoundBox& bb,
+        const bool removeInside,
+        List<surfaceFeatures::edgeStatus>& edgeStat
+    );
+
+    bool onLine(const point& p, const linePointRef& line);
+
+    //- Deletes all edges inside/outside bounding box from set.
+    void deleteEdges
+    (
+        const triSurface& surf,
+        const plane& cutPlane,
+        List<surfaceFeatures::edgeStatus>& edgeStat
+    );
+
+    void processHit
+    (
+        scalar& internalCloseness,
+        scalar& externalCloseness,
+        const label fi,
+        const triSurface& surf,
+        const point& start,
+        const point& p,
+        const point& end,
+        const vector& normal,
+        const vectorField& normals,
+        const List<pointIndexHit>& hitInfo
+    );
+
+    void drawHitProblem
+    (
+        const label fi,
+        const triSurface& surf,
+        const point& start,
+        const point& p,
+        const point& end,
+        const List<pointIndexHit>& hitInfo
+    );
+
+    //- Unmark non-manifold edges if individual triangles are not features
+    void unmarkBaffles
+    (
+        const triSurface& surf,
+        const scalar includedAngle,
+        List<surfaceFeatures::edgeStatus>& edgeStat
+    );
+
+    //- Divide into multiple normal bins
+    //  - return REGION if != 2 normals
+    //  - return REGION if 2 normals that make feature angle
+    //  - otherwise return NONE and set normals,bins
+    surfaceFeatures::edgeStatus checkFlatRegionEdge
+    (
+        const triSurface& surf,
+        const scalar tol,
+        const scalar includedAngle,
+        const label edgeI
+    );
+
+    void extractCloseness
+    (
+        const fileName &sFeatFileName,
+        const Time& runTime,
+        const triSurface &surf,
+        const bool writeVTK
+    );
+
+    void extractPointCloseness
+    (
+        const fileName &sFeatFileName,
+        const Time& runTime,
+        const triSurface &surf,
+        const bool writeVTK
+    );
+
+    void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os);
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C
new file mode 100644
index 000000000..e8660af3d
--- /dev/null
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C
@@ -0,0 +1,938 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    surfaceFeatureExtract
+
+Description
+    Extracts and writes surface features to file. All but the basic feature
+    extraction is WIP.
+
+    Curvature calculation is an implementation of the algorithm from:
+
+      "Estimating Curvatures and their Derivatives on Triangle Meshes"
+      by S. Rusinkiewicz
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeatureExtract.H"
+#include "Time.H"
+#include "meshTools.H"
+#include "tensor2D.H"
+#include "symmTensor2D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const Foam::scalar Foam::internalAngleTolerance(45);
+const Foam::scalar Foam::internalToleranceCosAngle
+(
+    cos(degToRad(180 - internalAngleTolerance))
+);
+
+const Foam::scalar Foam::externalAngleTolerance(10);
+const Foam::scalar Foam::externalToleranceCosAngle
+(
+    cos(degToRad(180 - externalAngleTolerance))
+);
+
+
+Foam::scalar Foam::calcVertexNormalWeight
+(
+    const triFace& f,
+    const label pI,
+    const vector& fN,
+    const pointField& points
+)
+{
+    label index = findIndex(f, pI);
+
+    if (index == -1)
+    {
+        FatalErrorInFunction
+            << "Point not in face" << abort(FatalError);
+    }
+
+    const vector e1 = points[f[index]] - points[f[f.fcIndex(index)]];
+    const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
+
+    return mag(fN)/(magSqr(e1)*magSqr(e2) + VSMALL);
+}
+
+
+Foam::point Foam::randomPointInPlane(const plane& p)
+{
+    // Perturb base point
+    const point& refPt = p.refPoint();
+
+    // ax + by + cz + d = 0
+    const FixedList<scalar, 4>& planeCoeffs = p.planeCoeffs();
+
+    const scalar perturbX = refPt.x() + 1e-3;
+    const scalar perturbY = refPt.y() + 1e-3;
+    const scalar perturbZ = refPt.z() + 1e-3;
+
+    if (mag(planeCoeffs[2]) < SMALL)
+    {
+        if (mag(planeCoeffs[1]) < SMALL)
+        {
+            const scalar x =
+                -1.0
+                *(
+                     planeCoeffs[3]
+                   + planeCoeffs[1]*perturbY
+                   + planeCoeffs[2]*perturbZ
+                 )/planeCoeffs[0];
+
+            return point
+            (
+                x,
+                perturbY,
+                perturbZ
+            );
+        }
+
+        const scalar y =
+            -1.0
+            *(
+                 planeCoeffs[3]
+               + planeCoeffs[0]*perturbX
+               + planeCoeffs[2]*perturbZ
+             )/planeCoeffs[1];
+
+        return point
+        (
+            perturbX,
+            y,
+            perturbZ
+        );
+    }
+    else
+    {
+        const scalar z =
+            -1.0
+            *(
+                 planeCoeffs[3]
+               + planeCoeffs[0]*perturbX
+               + planeCoeffs[1]*perturbY
+             )/planeCoeffs[2];
+
+        return point
+        (
+            perturbX,
+            perturbY,
+            z
+        );
+    }
+}
+
+
+Foam::triadField Foam::calcVertexCoordSys
+(
+    const triSurface& surf,
+    const vectorField& pointNormals
+)
+{
+    const pointField& points = surf.points();
+    const Map<label>& meshPointMap = surf.meshPointMap();
+
+    triadField pointCoordSys(points.size());
+
+    forAll(points, pI)
+    {
+        const point& pt = points[pI];
+        const vector& normal = pointNormals[meshPointMap[pI]];
+
+        if (mag(normal) < SMALL)
+        {
+            pointCoordSys[meshPointMap[pI]] = triad::unset;
+            continue;
+        }
+
+        plane p(pt, normal);
+
+        // Pick random point in plane
+        vector dir1 = pt - randomPointInPlane(p);
+        dir1 /= mag(dir1);
+
+        vector dir2 = dir1 ^ normal;
+        dir2 /= mag(dir2);
+
+        pointCoordSys[meshPointMap[pI]] = triad(dir1, dir2, normal);
+    }
+
+    return pointCoordSys;
+}
+
+
+Foam::vectorField Foam::calcVertexNormals(const triSurface& surf)
+{
+    // Weighted average of normals of faces attached to the vertex
+    // Weight = fA / (mag(e0)^2 * mag(e1)^2);
+
+    Info<< "Calculating vertex normals" << endl;
+
+    vectorField pointNormals(surf.nPoints(), Zero);
+
+    const pointField& points = surf.points();
+    const labelListList& pointFaces = surf.pointFaces();
+    const labelList& meshPoints = surf.meshPoints();
+
+    forAll(pointFaces, pI)
+    {
+        const labelList& pFaces = pointFaces[pI];
+
+        forAll(pFaces, fi)
+        {
+            const label facei = pFaces[fi];
+            const triFace& f = surf[facei];
+
+            vector fN = f.normal(points);
+
+            scalar weight = calcVertexNormalWeight
+            (
+                f,
+                meshPoints[pI],
+                fN,
+                points
+            );
+
+            pointNormals[pI] += weight*fN;
+        }
+
+        pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL;
+    }
+
+    return pointNormals;
+}
+
+
+Foam::triSurfacePointScalarField Foam::calcCurvature
+(
+    const word& name,
+    const Time& runTime,
+    const triSurface& surf,
+    const vectorField& pointNormals,
+    const triadField& pointCoordSys
+)
+{
+    Info<< "Calculating face curvature" << endl;
+
+    const pointField& points = surf.points();
+    const labelList& meshPoints = surf.meshPoints();
+    const Map<label>& meshPointMap = surf.meshPointMap();
+
+    triSurfacePointScalarField curvaturePointField
+    (
+        IOobject
+        (
+            name + ".curvature",
+            runTime.constant(),
+            "triSurface",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        scalarField(points.size(), 0.0)
+    );
+
+    List<symmTensor2D> pointFundamentalTensors
+    (
+        points.size(),
+        symmTensor2D::zero
+    );
+
+    scalarList accumulatedWeights(points.size(), 0.0);
+
+    forAll(surf, fi)
+    {
+        const triFace& f = surf[fi];
+        const edgeList fEdges = f.edges();
+
+        // Calculate the edge vectors and the normal differences
+        vectorField edgeVectors(f.size(), Zero);
+        vectorField normalDifferences(f.size(), Zero);
+
+        forAll(fEdges, feI)
+        {
+            const edge& e = fEdges[feI];
+
+            edgeVectors[feI] = e.vec(points);
+            normalDifferences[feI] =
+               pointNormals[meshPointMap[e[0]]]
+             - pointNormals[meshPointMap[e[1]]];
+        }
+
+        // Set up a local coordinate system for the face
+        const vector& e0 = edgeVectors[0];
+        const vector eN = f.normal(points);
+        const vector e1 = (e0 ^ eN);
+
+        if (magSqr(eN) < ROOTVSMALL)
+        {
+            continue;
+        }
+
+        triad faceCoordSys(e0, e1, eN);
+        faceCoordSys.normalize();
+
+        // Construct the matrix to solve
+        scalarSymmetricSquareMatrix T(3, 0);
+        scalarDiagonalMatrix Z(3, 0);
+
+        // Least Squares
+        for (label i = 0; i < 3; ++i)
+        {
+            scalar x = edgeVectors[i] & faceCoordSys[0];
+            scalar y = edgeVectors[i] & faceCoordSys[1];
+
+            T(0, 0) += sqr(x);
+            T(1, 0) += x*y;
+            T(1, 1) += sqr(x) + sqr(y);
+            T(2, 1) += x*y;
+            T(2, 2) += sqr(y);
+
+            scalar dndx = normalDifferences[i] & faceCoordSys[0];
+            scalar dndy = normalDifferences[i] & faceCoordSys[1];
+
+            Z[0] += dndx*x;
+            Z[1] += dndx*y + dndy*x;
+            Z[2] += dndy*y;
+        }
+
+        // Perform Cholesky decomposition and back substitution.
+        // Decomposed matrix is in T and solution is in Z.
+        LUsolve(T, Z);
+        symmTensor2D secondFundamentalTensor(Z[0], Z[1], Z[2]);
+
+        // Loop over the face points adding the contribution of the face
+        // curvature to the points.
+        forAll(f, fpI)
+        {
+            const label patchPointIndex = meshPointMap[f[fpI]];
+
+            const triad& ptCoordSys = pointCoordSys[patchPointIndex];
+
+            if (!ptCoordSys.set())
+            {
+                continue;
+            }
+
+            // Rotate faceCoordSys to ptCoordSys
+            tensor rotTensor = rotationTensor(ptCoordSys[2], faceCoordSys[2]);
+            triad rotatedFaceCoordSys = rotTensor & tensor(faceCoordSys);
+
+            // Project the face curvature onto the point plane
+
+            vector2D cmp1
+            (
+                ptCoordSys[0] & rotatedFaceCoordSys[0],
+                ptCoordSys[0] & rotatedFaceCoordSys[1]
+            );
+            vector2D cmp2
+            (
+                ptCoordSys[1] & rotatedFaceCoordSys[0],
+                ptCoordSys[1] & rotatedFaceCoordSys[1]
+            );
+
+            tensor2D projTensor
+            (
+                cmp1,
+                cmp2
+            );
+
+            symmTensor2D projectedFundamentalTensor
+            (
+                projTensor.x() & (secondFundamentalTensor & projTensor.x()),
+                projTensor.x() & (secondFundamentalTensor & projTensor.y()),
+                projTensor.y() & (secondFundamentalTensor & projTensor.y())
+            );
+
+            // Calculate weight
+            // TODO: Voronoi area weighting
+            scalar weight = calcVertexNormalWeight
+            (
+                f,
+                meshPoints[patchPointIndex],
+                f.normal(points),
+                points
+            );
+
+            // Sum contribution of face to this point
+            pointFundamentalTensors[patchPointIndex] +=
+                weight*projectedFundamentalTensor;
+
+            accumulatedWeights[patchPointIndex] += weight;
+        }
+
+        if (false)
+        {
+            Info<< "Points = " << points[f[0]] << " "
+                << points[f[1]] << " "
+                << points[f[2]] << endl;
+            Info<< "edgeVecs = " << edgeVectors[0] << " "
+                << edgeVectors[1] << " "
+                << edgeVectors[2] << endl;
+            Info<< "normDiff = " << normalDifferences[0] << " "
+                << normalDifferences[1] << " "
+                << normalDifferences[2] << endl;
+            Info<< "faceCoordSys = " << faceCoordSys << endl;
+            Info<< "T = " << T << endl;
+            Info<< "Z = " << Z << endl;
+        }
+    }
+
+    forAll(curvaturePointField, pI)
+    {
+        pointFundamentalTensors[pI] /= (accumulatedWeights[pI] + SMALL);
+
+        vector2D principalCurvatures = eigenValues(pointFundamentalTensors[pI]);
+
+        //scalar curvature =
+        //    (principalCurvatures[0] + principalCurvatures[1])/2;
+        scalar curvature = max
+        (
+            mag(principalCurvatures[0]),
+            mag(principalCurvatures[1])
+        );
+        //scalar curvature = principalCurvatures[0]*principalCurvatures[1];
+
+        curvaturePointField[meshPoints[pI]] = curvature;
+    }
+
+    return curvaturePointField;
+}
+
+
+bool Foam::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;
+}
+
+
+Foam::scalar Foam::calcProximityOfFeaturePoints
+(
+    const List<pointIndexHit>& hitList,
+    const scalar defaultCellSize
+)
+{
+    scalar minDist = defaultCellSize;
+
+    for
+    (
+        label hI1 = 0;
+        hI1 < hitList.size() - 1;
+        ++hI1
+    )
+    {
+        const pointIndexHit& pHit1 = hitList[hI1];
+
+        if (pHit1.hit())
+        {
+            for
+            (
+                label hI2 = hI1 + 1;
+                hI2 < hitList.size();
+                ++hI2
+            )
+            {
+                const pointIndexHit& pHit2 = hitList[hI2];
+
+                if (pHit2.hit())
+                {
+                    scalar curDist = mag(pHit1.hitPoint() - pHit2.hitPoint());
+
+                    minDist = min(curDist, minDist);
+                }
+            }
+        }
+    }
+
+    return minDist;
+}
+
+
+Foam::scalar Foam::calcProximityOfFeatureEdges
+(
+    const extendedFeatureEdgeMesh& efem,
+    const List<pointIndexHit>& hitList,
+    const scalar defaultCellSize
+)
+{
+    scalar minDist = defaultCellSize;
+
+    for
+    (
+        label hI1 = 0;
+        hI1 < hitList.size() - 1;
+        ++hI1
+    )
+    {
+        const pointIndexHit& pHit1 = hitList[hI1];
+
+        if (pHit1.hit())
+        {
+            const edge& e1 = efem.edges()[pHit1.index()];
+
+            for
+            (
+                label hI2 = hI1 + 1;
+                hI2 < hitList.size();
+                ++hI2
+            )
+            {
+                const pointIndexHit& pHit2 = hitList[hI2];
+
+                if (pHit2.hit())
+                {
+                    const edge& e2 = efem.edges()[pHit2.index()];
+
+                    // Don't refine if the edges are connected to each other
+                    if (!edgesConnected(e1, e2))
+                    {
+                        scalar curDist =
+                            mag(pHit1.hitPoint() - pHit2.hitPoint());
+
+                        minDist = min(curDist, minDist);
+                    }
+                }
+            }
+        }
+    }
+
+    return minDist;
+}
+
+
+void Foam::dumpBox(const treeBoundBox& bb, const fileName& fName)
+{
+    OFstream str(fName);
+
+    Info<< "Dumping bounding box " << bb << " as lines to obj file "
+        << str.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;
+    }
+}
+
+
+void Foam::deleteBox
+(
+    const triSurface& surf,
+    const treeBoundBox& bb,
+    const bool removeInside,
+    List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+    forAll(edgeStat, edgeI)
+    {
+        const point eMid = surf.edges()[edgeI].centre(surf.localPoints());
+
+        if (removeInside ? bb.contains(eMid) : !bb.contains(eMid))
+        {
+            edgeStat[edgeI] = surfaceFeatures::NONE;
+        }
+    }
+}
+
+
+bool Foam::onLine(const point& p, const linePointRef& line)
+{
+    const point& a = line.start();
+    const point& b = line.end();
+
+    if
+    (
+        ( p.x() < min(a.x(), b.x()) || p.x() > max(a.x(), b.x()) )
+     || ( p.y() < min(a.y(), b.y()) || p.y() > max(a.y(), b.y()) )
+     || ( p.z() < min(a.z(), b.z()) || p.z() > max(a.z(), b.z()) )
+    )
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
+void Foam::deleteEdges
+(
+    const triSurface& surf,
+    const plane& cutPlane,
+    List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+    const pointField& points = surf.points();
+    const labelList& meshPoints = surf.meshPoints();
+
+    forAll(edgeStat, edgeI)
+    {
+        const edge& e = surf.edges()[edgeI];
+        const point& p0 = points[meshPoints[e.start()]];
+        const point& p1 = points[meshPoints[e.end()]];
+        const linePointRef line(p0, p1);
+
+        // If edge does not intersect the plane, delete.
+        scalar intersect = cutPlane.lineIntersect(line);
+
+        point featPoint = intersect * (p1 - p0) + p0;
+
+        if (!onLine(featPoint, line))
+        {
+            edgeStat[edgeI] = surfaceFeatures::NONE;
+        }
+    }
+}
+
+
+void Foam::drawHitProblem
+(
+    const label fi,
+    const triSurface& surf,
+    const point& start,
+    const point& p,
+    const point& end,
+    const List<pointIndexHit>& hitInfo
+)
+{
+    Info<< nl << "# findLineAll did not hit its own face."
+        << nl << "# fi " << fi
+        << nl << "# start " << start
+        << nl << "# point " << p
+        << nl << "# end " << end
+        << nl << "# hitInfo " << hitInfo
+        << endl;
+
+    meshTools::writeOBJ(Info, start);
+    meshTools::writeOBJ(Info, p);
+    meshTools::writeOBJ(Info, end);
+
+    Info<< "l 1 2 3" << endl;
+
+    meshTools::writeOBJ(Info, surf.points()[surf[fi][0]]);
+    meshTools::writeOBJ(Info, surf.points()[surf[fi][1]]);
+    meshTools::writeOBJ(Info, surf.points()[surf[fi][2]]);
+
+    Info<< "f 4 5 6" << endl;
+
+    forAll(hitInfo, hI)
+    {
+        label hFI = hitInfo[hI].index();
+
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
+
+        Info<< "f "
+            << 3*hI + 7 << " "
+            << 3*hI + 8 << " "
+            << 3*hI + 9
+            << endl;
+    }
+}
+
+
+void Foam::unmarkBaffles
+(
+    const triSurface& surf,
+    const scalar includedAngle,
+    List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+    scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
+
+    const labelListList& edgeFaces = surf.edgeFaces();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() > 2)
+        {
+            label i0 = eFaces[0];
+            //const labelledTri& f0 = surf[i0];
+            const Foam::vector& n0 = surf.faceNormals()[i0];
+
+            //Pout<< "edge:" << edgeI << " n0:" << n0 << endl;
+
+            bool same = true;
+
+            for (label i = 1; i < eFaces.size(); i++)
+            {
+                //const labelledTri& f = surf[i];
+                const Foam::vector& n = surf.faceNormals()[eFaces[i]];
+
+                //Pout<< "    mag(n&n0): " << mag(n&n0) << endl;
+
+                if (mag(n&n0) < minCos)
+                {
+                    same = false;
+                    break;
+                }
+            }
+
+            if (same)
+            {
+                edgeStat[edgeI] = surfaceFeatures::NONE;
+            }
+        }
+    }
+}
+
+
+Foam::surfaceFeatures::edgeStatus Foam::checkFlatRegionEdge
+(
+    const triSurface& surf,
+    const scalar tol,
+    const scalar includedAngle,
+    const label edgeI
+)
+{
+    const edge& e = surf.edges()[edgeI];
+    const labelList& eFaces = surf.edgeFaces()[edgeI];
+
+    // Bin according to normal
+
+    DynamicList<Foam::vector> normals(2);
+    DynamicList<labelList> bins(2);
+
+    forAll(eFaces, eFacei)
+    {
+        const Foam::vector& n = surf.faceNormals()[eFaces[eFacei]];
+
+        // Find the normal in normals
+        label index = -1;
+        forAll(normals, normalI)
+        {
+            if (mag(n&normals[normalI]) > (1-tol))
+            {
+                index = normalI;
+                break;
+            }
+        }
+
+        if (index != -1)
+        {
+            bins[index].append(eFacei);
+        }
+        else if (normals.size() >= 2)
+        {
+            // Would be third normal. Mark as feature.
+            //Pout<< "** at edge:" << surf.localPoints()[e[0]]
+            //    << surf.localPoints()[e[1]]
+            //    << " have normals:" << normals
+            //    << " and " << n << endl;
+            return surfaceFeatures::REGION;
+        }
+        else
+        {
+            normals.append(n);
+            bins.append(labelList(1, eFacei));
+        }
+    }
+
+
+    // Check resulting number of bins
+    if (bins.size() == 1)
+    {
+        // Note: should check here whether they are two sets of faces
+        // that are planar or indeed 4 faces al coming together at an edge.
+        //Pout<< "** at edge:"
+        //    << surf.localPoints()[e[0]]
+        //    << surf.localPoints()[e[1]]
+        //    << " have single normal:" << normals[0]
+        //    << endl;
+        return surfaceFeatures::NONE;
+    }
+    else
+    {
+        // Two bins. Check if normals make an angle
+
+        //Pout<< "** at edge:"
+        //    << surf.localPoints()[e[0]]
+        //    << surf.localPoints()[e[1]] << nl
+        //    << "    normals:" << normals << nl
+        //    << "    bins   :" << bins << nl
+        //    << endl;
+
+        if (includedAngle >= 0)
+        {
+            scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
+
+            forAll(eFaces, i)
+            {
+                const Foam::vector& ni = surf.faceNormals()[eFaces[i]];
+                for (label j=i+1; j<eFaces.size(); j++)
+                {
+                    const Foam::vector& nj = surf.faceNormals()[eFaces[j]];
+                    if (mag(ni & nj) < minCos)
+                    {
+                        //Pout<< "have sharp feature between normal:" << ni
+                        //    << " and " << nj << endl;
+
+                        // Is feature. Keep as region or convert to
+                        // feature angle? For now keep as region.
+                        return surfaceFeatures::REGION;
+                    }
+                }
+            }
+        }
+
+
+        // So now we have two normals bins but need to make sure both
+        // bins have the same regions in it.
+
+         // 1. store + or - region number depending
+        //    on orientation of triangle in bins[0]
+        const labelList& bin0 = bins[0];
+        labelList regionAndNormal(bin0.size());
+        forAll(bin0, i)
+        {
+            const labelledTri& t = surf.localFaces()[eFaces[bin0[i]]];
+            int dir = t.edgeDirection(e);
+
+            if (dir > 0)
+            {
+                regionAndNormal[i] = t.region()+1;
+            }
+            else if (dir == 0)
+            {
+                FatalErrorInFunction
+                    << exit(FatalError);
+            }
+            else
+            {
+                regionAndNormal[i] = -(t.region()+1);
+            }
+        }
+
+        // 2. check against bin1
+        const labelList& bin1 = bins[1];
+        labelList regionAndNormal1(bin1.size());
+        forAll(bin1, i)
+        {
+            const labelledTri& t = surf.localFaces()[eFaces[bin1[i]]];
+            int dir = t.edgeDirection(e);
+
+            label myRegionAndNormal;
+            if (dir > 0)
+            {
+                myRegionAndNormal = t.region()+1;
+            }
+            else
+            {
+                myRegionAndNormal = -(t.region()+1);
+            }
+
+            regionAndNormal1[i] = myRegionAndNormal;
+
+            label index = findIndex(regionAndNormal, -myRegionAndNormal);
+            if (index == -1)
+            {
+                // Not found.
+                //Pout<< "cannot find region " << myRegionAndNormal
+                //    << " in regions " << regionAndNormal << endl;
+
+                return surfaceFeatures::REGION;
+            }
+        }
+
+        // Passed all checks, two normal bins with the same contents.
+        //Pout<< "regionAndNormal:" << regionAndNormal << endl;
+        //Pout<< "myRegionAndNormal:" << regionAndNormal1 << endl;
+
+        return surfaceFeatures::NONE;
+    }
+}
+
+
+void Foam::extractCloseness
+(
+    const fileName &sFeatFileName,
+    const Time& runTime,
+    const triSurface &surf,
+    const bool writeVTK
+);
+
+
+void Foam::extractPointCloseness
+(
+    const fileName &sFeatFileName,
+    const Time& runTime,
+    const triSurface &surf,
+    const bool writeVTK
+);
+
+
+void Foam::writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
+{
+    os  << "    points : " << fem.points().size() << nl
+        << "    of which" << nl
+        << "        convex             : "
+        << fem.concaveStart() << nl
+        << "        concave            : "
+        << (fem.mixedStart()-fem.concaveStart()) << nl
+        << "        mixed              : "
+        << (fem.nonFeatureStart()-fem.mixedStart()) << nl
+        << "        non-feature        : "
+        << (fem.points().size()-fem.nonFeatureStart()) << nl
+        << "    edges  : " << fem.edges().size() << nl
+        << "    of which" << nl
+        << "        external edges     : "
+        << fem.internalStart() << nl
+        << "        internal edges     : "
+        << (fem.flatStart()- fem.internalStart()) << nl
+        << "        flat edges         : "
+        << (fem.openStart()- fem.flatStart()) << nl
+        << "        open edges         : "
+        << (fem.multipleStart()- fem.openStart()) << nl
+        << "        multiply connected : "
+        << (fem.edges().size()- fem.multipleStart()) << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
index 793d32fc7..cd57e5dec 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,9 +66,9 @@ bool Foam::triSurfaceSearch::checkUniqueHit
 
             if (pointFacei != currHit.index())
             {
-                forAll(hits, hI)
+                forAll(hits, hi)
                 {
-                    const pointIndexHit& hit = hits[hI];
+                    const pointIndexHit& hit = hits[hi];
 
                     if (hit.index() == pointFacei)
                     {
@@ -89,15 +89,15 @@ bool Foam::triSurfaceSearch::checkUniqueHit
 
         const labelList& edgeFaces = surface().edgeFaces()[edgeI];
 
-        forAll(edgeFaces, fI)
+        forAll(edgeFaces, fi)
         {
-            const label edgeFacei = edgeFaces[fI];
+            const label edgeFacei = edgeFaces[fi];
 
             if (edgeFacei != currHit.index())
             {
-                forAll(hits, hI)
+                forAll(hits, hi)
                 {
-                    const pointIndexHit& hit = hits[hI];
+                    const pointIndexHit& hit = hits[hi];
 
                     if (hit.index() == edgeFacei)
                     {
@@ -294,7 +294,7 @@ void Foam::triSurfaceSearch::findNearest
 
     forAll(samples, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findNearest
+        info[i] = octree.findNearest
         (
             samples[i],
             nearestDistSqr[i],
@@ -335,7 +335,7 @@ void Foam::triSurfaceSearch::findLine
 
     forAll(start, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findLine
+        info[i] = octree.findLine
         (
             start[i],
             end[i]
@@ -362,7 +362,7 @@ void Foam::triSurfaceSearch::findLineAny
 
     forAll(start, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
+        info[i] = octree.findLineAny
         (
             start[i],
             end[i]
@@ -394,7 +394,7 @@ void Foam::triSurfaceSearch::findLineAll
 
     treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
 
-    forAll(start, pointi)
+    forAll(start, i)
     {
         hits.clear();
         shapeMask.clear();
@@ -404,25 +404,17 @@ void Foam::triSurfaceSearch::findLineAll
             // See if any intersection between pt and end
             pointIndexHit inter = octree.findLine
             (
-                start[pointi],
-                end[pointi],
+                start[i],
+                end[i],
                 allIntersectOp
             );
 
             if (inter.hit())
             {
-                vector lineVec = end[pointi] - start[pointi];
+                vector lineVec = end[i] - start[i];
                 lineVec /= mag(lineVec) + VSMALL;
 
-                if
-                (
-                    checkUniqueHit
-                    (
-                        inter,
-                        hits,
-                        lineVec
-                    )
-                )
+                if (checkUniqueHit(inter, hits, lineVec))
                 {
                     hits.append(inter);
                 }
@@ -435,7 +427,7 @@ void Foam::triSurfaceSearch::findLineAll
             }
         }
 
-        info[pointi].transfer(hits);
+        info[i].transfer(hits);
     }
 
     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
diff --git a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H
index e6a6a91c0..1a3f07d21 100644
--- a/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H
+++ b/src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -170,7 +170,6 @@ public:
             const bool isNodeValues,
             const bool verbose = false
         ) const;
-
 };
 
 
diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C
index 51ee83a03..fc9bf1418 100644
--- a/src/triSurface/triSurface/triSurface.C
+++ b/src/triSurface/triSurface/triSurface.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1031,6 +1031,19 @@ Foam::triSurface Foam::triSurface::subsetMesh
 }
 
 
+Foam::faceList Foam::triSurface::faces() const
+{
+    faceList faces(size());
+
+    forAll(*this, facei)
+    {
+        faces[facei] = operator[](facei).triFaceFace();
+    }
+
+    return faces;
+}
+
+
 void Foam::triSurface::write
 (
     const fileName& name,
diff --git a/src/triSurface/triSurface/triSurface.H b/src/triSurface/triSurface/triSurface.H
index 2e9ce52d1..a6082c385 100644
--- a/src/triSurface/triSurface/triSurface.H
+++ b/src/triSurface/triSurface/triSurface.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -393,6 +393,12 @@ public:
             ) const;
 
 
+        // Conversion
+
+            //- Return the list of triangles as a faceList
+            faceList faces() const;
+
+
         // Write
 
             //- Write to Ostream in simple FOAM format