diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 3ae006e69d..1ea773f16b 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -31,817 +31,31 @@ Description
     Extracts and writes surface features to file. All but the basic feature
     extraction is a work-in-progress.
 
-    The curvature calculation is an implementation of the algorithm from:
-    \verbatim
-        "Estimating Curvatures and their Derivatives on Triangle Meshes"
-        by S. Rusinkiewicz
-        3DPVT'04 Proceedings of the 3D Data Processing,
-        Visualization, and Transmission, 2nd International Symposium
-        Pages 486-493
-        http://gfx.cs.princeton.edu/pubs/_2004_ECA/curvpaper.pdf
-    \endverbatim
-
 \*---------------------------------------------------------------------------*/
 
 #include "argList.H"
 #include "Time.H"
 #include "triSurface.H"
+#include "triSurfaceTools.H"
+#include "edgeMeshTools.H"
 #include "surfaceFeaturesExtraction.H"
 #include "surfaceIntersection.H"
 #include "featureEdgeMesh.H"
 #include "extendedFeatureEdgeMesh.H"
 #include "treeBoundBox.H"
 #include "meshTools.H"
-#include "OFstream.H"
+#include "OBJstream.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 "triSurfaceLoader.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);
-}
-
-
-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 arbitrary point in plane
-        vector dir1 = pt - p.somePointInPlane(1e-3);
-        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;
-}
-
-
-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 (!e1.connects(e2))
-                    {
-                        scalar curDist =
-                            mag(pHit1.hitPoint() - pHit2.hitPoint());
-
-                        minDist = min(curDist, minDist);
-                    }
-                }
-            }
-        }
-    }
-
-    return minDist;
-}
-
-
-void dumpBox(const treeBoundBox& bb, const fileName& fName)
-{
-    OFstream os(fName);
-
-    Info<< "Dumping bounding box " << bb << " as lines to obj file "
-        << os.name() << endl;
-
-    meshTools::writeOBJ(os, bb);
-}
-
-
-// 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
@@ -874,7 +88,6 @@ int main(int argc, char *argv[])
         {
             continue;
         }
-
         const dictionary& surfaceDict = iter().dict();
 
         if (!surfaceDict.found("extractionMethod"))
@@ -965,7 +178,7 @@ int main(int argc, char *argv[])
         surf.writeStats(Info);
         Info<< nl;
 
-        // need plain faces if outputting VTK format
+        // Need a copy as plain faces if outputting VTK format
         faceList faces;
         if (writeVTK)
         {
@@ -976,15 +189,15 @@ int main(int argc, char *argv[])
             }
         }
 
-
-        // Either construct features from surface & featureAngle or read set.
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
+        //
+        // Extract features using the preferred extraction method
+        //
         autoPtr<surfaceFeatures> features = extractor().features(surf);
 
         // Trim set
         // ~~~~~~~~
 
+        // Option: "trimFeatures" (dictionary)
         if (surfaceDict.isDict("trimFeatures"))
         {
             const dictionary& trimDict = surfaceDict.subDict("trimFeatures");
@@ -1015,7 +228,6 @@ int main(int argc, char *argv[])
             }
         }
 
-
         // Subset
         // ~~~~~~
 
@@ -1035,51 +247,51 @@ int main(int argc, char *argv[])
             {
                 treeBoundBox bb(subsetDict.lookup("insideBox")());
 
-                Info<< "Removing all edges outside bb " << bb << endl;
-                dumpBox(bb, "subsetBox.obj");
+                Info<< "Subset edges inside box " << bb << endl;
+                features().subsetBox(edgeStat, bb);
 
-                deleteBox(surf, bb, false, edgeStat);
+                {
+                    OBJstream os("subsetBox.obj");
+
+                    Info<< "Dumping bounding box " << bb
+                        << " as lines to obj file "
+                        << os.name() << endl;
+
+                    os.write(bb);
+                }
             }
             // Suboption: "outsideBox"
             else if (subsetDict.found("outsideBox"))
             {
                 treeBoundBox bb(subsetDict.lookup("outsideBox")());
 
-                Info<< "Removing all edges inside bb " << bb << endl;
-                dumpBox(bb, "deleteBox.obj");
+                Info<< "Exclude edges outside box " << bb << endl;
+                features().excludeBox(edgeStat, bb);
 
-                deleteBox(surf, bb, true, edgeStat);
+                {
+                    OBJstream os("deleteBox.obj");
+
+                    Info<< "Dumping bounding box " << bb
+                        << " as lines to obj file "
+                        << os.name() << endl;
+
+                    os.write(bb);
+                }
             }
 
-            const Switch nonManifoldEdges =
-                subsetDict.lookupOrDefault<Switch>("nonManifoldEdges", "yes");
-
-            if (!nonManifoldEdges)
+            // Suboption: "nonManifoldEdges" (false: remove non-manifold edges)
+            if (!subsetDict.lookupOrDefault<bool>("nonManifoldEdges", true))
             {
                 Info<< "Removing all non-manifold edges"
                     << " (edges with > 2 connected faces) unless they"
                     << " cross multiple regions" << endl;
 
-                forAll(edgeStat, edgeI)
-                {
-                    const labelList& eFaces = surf.edgeFaces()[edgeI];
-
-                    if
-                    (
-                        eFaces.size() > 2
-                     && edgeStat[edgeI] == surfaceFeatures::REGION
-                     && (eFaces.size() % 2) == 0
-                    )
-                    {
-                        edgeStat[edgeI] = checkFlatRegionEdge
-                        (
-                            surf,
-                            1e-5,   //tol,
-                            extractor().includedAngle(),
-                            edgeI
-                        );
-                    }
-                }
+                features().checkFlatRegionEdge
+                (
+                    edgeStat,
+                    1e-5,   // tol
+                    extractor().includedAngle()
+                );
             }
 
             // Suboption: "openEdges" (false: remove open edges)
@@ -1088,13 +300,7 @@ int main(int argc, char *argv[])
                 Info<< "Removing all open edges"
                     << " (edges with 1 connected face)" << endl;
 
-                forAll(edgeStat, edgeI)
-                {
-                    if (surf.edgeFaces()[edgeI].size() == 1)
-                    {
-                        edgeStat[edgeI] = surfaceFeatures::NONE;
-                    }
-                }
+                features().excludeOpen(edgeStat);
             }
 
             // Suboption: "plane"
@@ -1102,47 +308,49 @@ int main(int argc, char *argv[])
             {
                 plane cutPlane(subsetDict.lookup("plane")());
 
-                deleteEdges(surf, cutPlane, edgeStat);
+                Info<< "Only include feature edges that intersect the plane"
+                    << " with normal " << cutPlane.normal()
+                    << " and base point " << cutPlane.refPoint() << endl;
 
-                Info<< "Only edges that intersect the plane with normal "
-                    << cutPlane.normal()
-                    << " and base point " << cutPlane.refPoint()
-                    << " will be included as feature edges."<< endl;
+                features().subsetPlane(edgeStat, cutPlane);
             }
         }
 
-
         surfaceFeatures newSet(surf);
         newSet.setFromStatus(edgeStat, extractor().includedAngle());
 
-        Info<< nl
-            << "Initial feature set:" << nl
-            << "    feature points : " << newSet.featurePoints().size() << nl
-            << "    feature edges  : " << newSet.featureEdges().size() << nl
-            << "    of which" << nl
-            << "        region edges   : " << newSet.nRegionEdges() << nl
-            << "        external edges : " << newSet.nExternalEdges() << nl
-            << "        internal edges : " << newSet.nInternalEdges() << nl
-            << endl;
-
-        //if (writeObj)
-        //{
-        //    newSet.writeObj("final");
-        //}
+        Info<< nl << "Initial ";
+        newSet.writeStats(Info);
 
         boolList surfBaffleRegions(surf.patches().size(), false);
-
-        wordList surfBaffleNames;
-        surfaceDict.readIfPresent("baffles", surfBaffleNames);
-
-        forAll(surf.patches(), pI)
+        if (surfaceDict.found("baffles"))
         {
-            const word& name = surf.patches()[pI].name();
+            wordReList baffleSelect(surfaceDict.lookup("baffles"));
 
-            if (findIndex(surfBaffleNames, name) != -1)
+            wordList patchNames(surf.patches().size());
+            forAll(surf.patches(), patchi)
             {
-                Info<< "Adding baffle region " << name << endl;
-                surfBaffleRegions[pI] = true;
+                patchNames[patchi] = surf.patches()[patchi].name();
+            }
+
+            labelList indices = findStrings(baffleSelect, patchNames);
+            forAll(indices, patchi)
+            {
+                surfBaffleRegions[patchi] = true;
+            }
+
+            if (indices.size())
+            {
+                Info<< "Adding " << indices.size() << " baffle regions: (";
+
+                forAll(surfBaffleRegions, patchi)
+                {
+                    if (surfBaffleRegions[patchi])
+                    {
+                        Info<< ' ' << patchNames[patchi];
+                    }
+                }
+                Info<< " )" << nl << nl;
             }
         }
 
@@ -1175,7 +383,7 @@ int main(int argc, char *argv[])
                 )
             );
             Info<< "Read " << addFeMesh.name() << nl;
-            writeStats(addFeMesh, Info);
+            edgeMeshTools::writeStats(Info, addFeMesh);
 
             feMesh.add(addFeMesh);
         }
@@ -1209,15 +417,14 @@ int main(int argc, char *argv[])
                 sizeInfo[0] = addMesh.points().size();
                 sizeInfo[1] = addMesh.edges().size();
             }
-            Info<< "Self intersection:" << nl
-                    << "    points : " << sizeInfo[0] << nl
-                    << "    edges  : " << sizeInfo[1] << nl;
+            Info<< nl
+                << "Self intersection:" << nl
+                << "    points : " << sizeInfo[0] << nl
+                << "    edges  : " << sizeInfo[1] << nl;
         }
 
-
-        Info<< nl
-            << "Final feature set:" << nl;
-        writeStats(feMesh, Info);
+        Info<< nl << "Final ";
+        edgeMeshTools::writeStats(Info, feMesh);
 
         Info<< nl << "Writing extendedFeatureEdgeMesh to "
             << feMesh.objectPath() << endl;
@@ -1259,253 +466,16 @@ int main(int argc, char *argv[])
         // Option: "closeness"
         if (surfaceDict.lookupOrDefault<bool>("closeness", false))
         {
-            Info<< nl << "Extracting internal and external closeness of "
-                << "surface." << endl;
-
-            triSurfaceMesh searchSurf
-            (
-                IOobject
+            Pair<tmp<scalarField>> tcloseness =
+                triSurfaceTools::writeCloseness
                 (
-                    outputName + ".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)
+                    outputName,
+                    surf,
+                    45,  // internalAngleTolerance
+                    10   // externalAngleTolerance
                 );
 
-            scalar internalAngleTolerance = 45;
-            scalar internalToleranceCosAngle =
-                Foam::cos
-                (
-                    degToRad(180 - internalAngleTolerance)
-                );
-
-            Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle
-                << nl
-                << "internalToleranceCosAngle: " << internalToleranceCosAngle
-                << endl;
-
-            // Info<< "span " << span << endl;
-
-            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];
-
-                if (hitInfo.size() < 1)
-                {
-                    drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
-
-                    // FatalErrorInFunction
-                    //     << "findLineAll did not hit its own face."
-                    //     << exit(FatalError);
-                }
-                else if (hitInfo.size() == 1)
-                {
-                    if (!hitInfo[0].hit())
-                    {
-                        // FatalErrorInFunction
-                        //     << "findLineAll did not hit any face."
-                        //     << exit(FatalError);
-                    }
-                    else if (hitInfo[0].index() != fI)
-                    {
-                        drawHitProblem
-                        (
-                            fI,
-                            surf,
-                            start,
-                            faceCentres,
-                            end,
-                            hitInfo
-                        );
-
-                        // FatalErrorInFunction
-                        //     << "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
-                        );
-
-                        // FatalErrorInFunction
-                        //     << "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
-                (
-                    outputName + ".internalCloseness",
-                    runTime.constant(),
-                    "triSurface",
-                    runTime,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                surf,
-                dimLength,
-                internalCloseness
-            );
-
-            internalClosenessField.write();
-
-            triSurfaceScalarField externalClosenessField
-            (
-                IOobject
-                (
-                    outputName + ".externalCloseness",
-                    runTime.constant(),
-                    "triSurface",
-                    runTime,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                surf,
-                dimLength,
-                externalCloseness
-            );
-
-            externalClosenessField.write();
-
             if (writeVTK)
             {
                 vtkSurfaceWriter().write
@@ -1518,7 +488,7 @@ int main(int argc, char *argv[])
                         faces
                     ),
                     "internalCloseness",                // fieldName
-                    internalCloseness,
+                    tcloseness[0](),
                     false,                              // isNodeValues
                     true                                // verbose
                 );
@@ -1533,7 +503,7 @@ int main(int argc, char *argv[])
                         faces
                     ),
                     "externalCloseness",                // fieldName
-                    externalCloseness,
+                    tcloseness[1](),
                     false,                              // isNodeValues
                     true                                // verbose
                 );
@@ -1543,22 +513,13 @@ int main(int argc, char *argv[])
         // Option: "curvature"
         if (surfaceDict.lookupOrDefault<bool>("curvature", false))
         {
-            Info<< nl << "Extracting curvature of surface at the points."
-                << endl;
-
-            vectorField pointNormals = calcVertexNormals(surf);
-            triadField pointCoordSys = calcVertexCoordSys(surf, pointNormals);
-
-            triSurfacePointScalarField k = calcCurvature
-            (
-                outputName,
-                runTime,
-                surf,
-                pointNormals,
-                pointCoordSys
-            );
-
-            k.write();
+            tmp<scalarField> tcurvatureField =
+                triSurfaceTools::writeCurvature
+                (
+                    runTime,
+                    outputName,
+                    surf
+                );
 
             if (writeVTK)
             {
@@ -1572,7 +533,7 @@ int main(int argc, char *argv[])
                         faces
                     ),
                     "curvature",                        // fieldName
-                    k,
+                    tcurvatureField(),
                     true,                               // isNodeValues
                     true                                // verbose
                 );
@@ -1582,64 +543,15 @@ int main(int argc, char *argv[])
         // Option: "featureProximity"
         if (surfaceDict.lookupOrDefault<bool>("featureProximity", false))
         {
-            Info<< nl << "Extracting proximity of close feature points and "
-                << "edges to the surface" << endl;
-
-            const scalar searchDistance =
-                readScalar(surfaceDict.lookup("maxFeatureProximity"));
-
-            scalarField featureProximity(surf.size(), searchDistance);
-
-            forAll(surf, fI)
-            {
-                const triPointRef& tri = surf[fI].tri(surf.points());
-                const point& triCentre = tri.circumCentre();
-
-                const scalar radiusSqr = min
+            tmp<scalarField> tproximity =
+                edgeMeshTools::writeFeatureProximity
                 (
-                    sqr(4*tri.circumRadius()),
-                    sqr(searchDistance)
-                );
-
-                List<pointIndexHit> hitList;
-
-                feMesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
-
-                featureProximity[fI] =
-                    calcProximityOfFeatureEdges
-                    (
-                        feMesh,
-                        hitList,
-                        featureProximity[fI]
-                    );
-
-                feMesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
-
-                featureProximity[fI] =
-                    calcProximityOfFeaturePoints
-                    (
-                        hitList,
-                        featureProximity[fI]
-                    );
-            }
-
-            triSurfaceScalarField featureProximityField
-            (
-                IOobject
-                (
-                    outputName + ".featureProximity",
-                    runTime.constant(),
-                    "triSurface",
                     runTime,
-                    IOobject::NO_READ,
-                    IOobject::NO_WRITE
-                ),
-                surf,
-                dimLength,
-                featureProximity
-            );
-
-            featureProximityField.write();
+                    outputName,
+                    feMesh,
+                    surf,
+                    readScalar(surfaceDict.lookup("maxFeatureProximity"))
+                );
 
             if (writeVTK)
             {
@@ -1653,7 +565,7 @@ int main(int argc, char *argv[])
                         faces
                     ),
                     "featureProximity",                 // fieldName
-                    featureProximity,
+                    tproximity(),
                     false,                              // isNodeValues
                     true                                // verbose
                 );
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index b77d0bdbc9..f9af11d9be 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -25,6 +25,12 @@ $(edgeMeshFormats)/starcd/STARCDedgeFormatRunTime.C
 $(edgeMeshFormats)/vtk/VTKedgeFormat.C
 $(edgeMeshFormats)/vtk/VTKedgeFormatRunTime.C
 
+
+edgeMeshTools = $(em)/edgeMeshTools
+$(edgeMeshTools)/edgeMeshTools.C
+$(edgeMeshTools)/edgeMeshFeatureProximity.C
+
+
 $(em)/featureEdgeMesh/featureEdgeMesh.C
 
 eem = $(em)/extendedEdgeMesh
@@ -215,6 +221,8 @@ triSurface/triangleFuncs/triangleFuncs.C
 triSurface/surfaceFeatures/surfaceFeatures.C
 triSurface/triSurfaceLoader/triSurfaceLoader.C
 triSurface/triSurfaceTools/triSurfaceTools.C
+triSurface/triSurfaceTools/triSurfaceCloseness.C
+triSurface/triSurfaceTools/triSurfaceCurvature.C
 triSurface/triSurfaceTools/geompack/geompack.C
 triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
 
diff --git a/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshFeatureProximity.C b/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshFeatureProximity.C
new file mode 100644
index 0000000000..c5e9e883e7
--- /dev/null
+++ b/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshFeatureProximity.C
@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "edgeMeshTools.H"
+
+#include "extendedEdgeMesh.H"
+#include "triSurface.H"
+#include "triSurfaceFields.H"
+#include "pointIndexHit.H"
+#include "MeshedSurface.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+static 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 edgeMesh& emesh,
+    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 = emesh.edges()[pHit1.index()];
+
+            for
+            (
+                label hI2 = hI1 + 1;
+                hI2 < hitList.size();
+                ++hI2
+            )
+            {
+                const pointIndexHit& pHit2 = hitList[hI2];
+
+                if (pHit2.hit())
+                {
+                    const edge& e2 = emesh.edges()[pHit2.index()];
+
+                    // Don't refine if the edges are connected to each other
+                    if (!e1.connects(e2))
+                    {
+                        scalar curDist =
+                            mag(pHit1.hitPoint() - pHit2.hitPoint());
+
+                        minDist = min(curDist, minDist);
+                    }
+                }
+            }
+        }
+    }
+
+    return minDist;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::edgeMeshTools::featureProximity
+(
+    const extendedEdgeMesh& emesh,
+    const triSurface& surf,
+    const scalar searchDistance
+)
+{
+    tmp<scalarField> tfld(new scalarField(surf.size(), searchDistance));
+    scalarField& featureProximity = tfld.ref();
+
+    Info<< "Extracting proximity of close feature points and "
+        << "edges to the surface" << endl;
+
+    forAll(surf, fI)
+    {
+        const triPointRef& tri = surf[fI].tri(surf.points());
+        const point& triCentre = tri.circumCentre();
+
+        const scalar radiusSqr = min
+        (
+            sqr(4*tri.circumRadius()),
+            sqr(searchDistance)
+        );
+
+        List<pointIndexHit> hitList;
+
+        emesh.allNearestFeatureEdges(triCentre, radiusSqr, hitList);
+
+        featureProximity[fI] =
+            calcProximityOfFeatureEdges
+            (
+                emesh,
+                hitList,
+                featureProximity[fI]
+            );
+
+        emesh.allNearestFeaturePoints(triCentre, radiusSqr, hitList);
+
+        featureProximity[fI] =
+            calcProximityOfFeaturePoints
+            (
+                hitList,
+                featureProximity[fI]
+            );
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::edgeMeshTools::writeFeatureProximity
+(
+    const Time& runTime,
+    const word& basename,
+    const extendedEdgeMesh& emesh,
+    const triSurface& surf,
+    const scalar searchDistance
+)
+{
+    Info<< nl << "Extracting curvature of surface at the points."
+        << endl;
+
+
+    tmp<scalarField> tfld =
+        edgeMeshTools::featureProximity(emesh, surf, searchDistance);
+    scalarField& featureProximity = tfld.ref();
+
+    triSurfaceScalarField outputField
+    (
+        IOobject
+        (
+            basename + ".featureProximity",
+            runTime.constant(),
+            "triSurface",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        scalarField()
+    );
+
+    outputField.swap(featureProximity);
+    outputField.write();
+    outputField.swap(featureProximity);
+
+    return tfld;
+}
+
+// ************************************************************************* //
diff --git a/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshTools.C b/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshTools.C
new file mode 100644
index 0000000000..eac159b8e9
--- /dev/null
+++ b/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshTools.C
@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "edgeMeshTools.H"
+
+#include "extendedFeatureEdgeMesh.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::edgeMeshTools::writeStats
+(
+    Ostream& os,
+    const extendedFeatureEdgeMesh& emesh
+)
+{
+    os  << "Feature set:" << nl
+        << "    points : " << emesh.points().size() << nl
+        << "    of which" << nl
+        << "        convex             : "
+        << emesh.concaveStart() << nl
+        << "        concave            : "
+        << (emesh.mixedStart()-emesh.concaveStart()) << nl
+        << "        mixed              : "
+        << (emesh.nonFeatureStart()-emesh.mixedStart()) << nl
+        << "        non-feature        : "
+        << (emesh.points().size()-emesh.nonFeatureStart()) << nl
+        << "    edges  : " << emesh.edges().size() << nl
+        << "    of which" << nl
+        << "        external edges     : "
+        << emesh.internalStart() << nl
+        << "        internal edges     : "
+        << (emesh.flatStart()- emesh.internalStart()) << nl
+        << "        flat edges         : "
+        << (emesh.openStart()- emesh.flatStart()) << nl
+        << "        open edges         : "
+        << (emesh.multipleStart()- emesh.openStart()) << nl
+        << "        multiply connected : "
+        << (emesh.edges().size()- emesh.multipleStart()) << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshTools.H b/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshTools.H
new file mode 100644
index 0000000000..aab9216f34
--- /dev/null
+++ b/src/meshTools/edgeMesh/edgeMeshTools/edgeMeshTools.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Namespace
+    Foam::edgeMeshTools
+
+Description
+    Collection of static functions related to edgeMesh features.
+
+SourceFiles
+    edgeMeshTools.C
+    edgeMeshFeatureProximity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef edgeMeshTools_H
+#define edgeMeshTools_H
+
+#include "tmp.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class edgeMesh;
+class extendedEdgeMesh;
+class extendedFeatureEdgeMesh;
+class triSurface;
+class Time;
+class Ostream;
+
+/*---------------------------------------------------------------------------*\
+                      Namespace edgeMeshTools Declaration
+\*---------------------------------------------------------------------------*/
+
+namespace edgeMeshTools
+{
+    //- Write some information
+    void writeStats
+    (
+        Ostream& os,
+        const extendedFeatureEdgeMesh& emesh
+    );
+
+
+    //- Calculate proximity of the features to the surface
+    tmp<scalarField> featureProximity
+    (
+        const extendedEdgeMesh& emesh,
+        const triSurface& surf,
+        const scalar searchDistance
+    );
+
+
+    //- Calculate proximity of the features to the surface and write the field
+    tmp<scalarField> writeFeatureProximity
+    (
+        const Time& runTime,
+        const word& basename,
+        const extendedEdgeMesh& emesh,
+        const triSurface& surf,
+        const scalar searchDistance
+    );
+
+} // End namespace edgeMeshTools
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
index 914287df13..71a4b81779 100644
--- a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
+++ b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H
@@ -94,23 +94,22 @@ public:
 
     enum pointStatus
     {
-        CONVEX,     // Fully convex point (w.r.t normals)
-        CONCAVE,    // Fully concave point
-        MIXED,      // A point surrounded by both convex and concave edges
-        NONFEATURE  // Not a feature point
+        CONVEX,     //!< Fully convex point (w.r.t normals)
+        CONCAVE,    //!< Fully concave point
+        MIXED,      //!< A point surrounded by both convex and concave edges
+        NONFEATURE  //!< Not a feature point
     };
 
     static const Foam::NamedEnum<pointStatus, 4> pointStatusNames_;
 
     enum edgeStatus
     {
-        EXTERNAL,   // "Convex" edge
-        INTERNAL,   // "Concave" edge
-        FLAT,       // Neither concave or convex, on a flat surface
-        OPEN,       // i.e. only connected to one face
-        MULTIPLE,   // Multiply connected (connected to more than two faces)
-        NONE        // Not a classified feature edge (consistency with
-                    // surfaceFeatures)
+        EXTERNAL,   //!< "Convex" edge
+        INTERNAL,   //!< "Concave" edge
+        FLAT,       //!< Neither concave or convex, on a flat surface
+        OPEN,       //!< Only connected to a single face
+        MULTIPLE,   //!< Multiply connected (connected to more than two faces)
+        NONE        //!< Unclassified (consistency with surfaceFeatures)
     };
 
     static const Foam::NamedEnum<edgeStatus, 6> edgeStatusNames_;
@@ -118,10 +117,10 @@ public:
     //- Normals point to the outside
     enum sideVolumeType
     {
-        INSIDE  = 0,  // mesh inside
-        OUTSIDE = 1,  // mesh outside
-        BOTH    = 2,  // e.g. a baffle
-        NEITHER = 3   // not sure when this may be used
+        INSIDE  = 0,  //!< mesh inside
+        OUTSIDE = 1,  //!< mesh outside
+        BOTH    = 2,  //!< e.g. a baffle
+        NEITHER = 3   //!< not sure when this may be used
     };
 
     static const Foam::NamedEnum<sideVolumeType, 4> sideVolumeTypeNames_;
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
index 055bf8fd63..c727121a60 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,7 +42,30 @@ namespace Foam
     defineTypeNameAndDebug(surfaceFeatures, 0);
 
     const scalar surfaceFeatures::parallelTolerance = sin(degToRad(1.0));
+
+
+//! \cond fileScope
+//  Check if the point is on the line
+static bool onLine(const Foam::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;
 }
+//! \endcond
+
+} // End namespace Foam
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -442,6 +465,181 @@ Foam::surfaceFeatures::labelScalar Foam::surfaceFeatures::walkSegment
 }
 
 
+//- 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
+//
+// This has been relocated from surfaceFeatureExtract and could be cleaned
+// up some more.
+//
+Foam::surfaceFeatures::edgeStatus
+Foam::surfaceFeatures::surfaceFeatures::checkFlatRegionEdge
+(
+    const scalar tol,
+    const scalar includedAngle,
+    const label edgeI
+) const
+{
+    const triSurface& surf = surf_;
+
+    const edge& e = surf.edges()[edgeI];
+    const labelList& eFaces = surf.edgeFaces()[edgeI];
+
+    // Bin according to normal
+
+    DynamicList<vector> normals(2);
+    DynamicList<labelList> bins(2);
+
+    forAll(eFaces, eFacei)
+    {
+        const 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 vector& ni = surf.faceNormals()[eFaces[i]];
+                for (label j=i+1; j<eFaces.size(); j++)
+                {
+                    const 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;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::surfaceFeatures::surfaceFeatures(const triSurface& surf)
@@ -829,9 +1027,127 @@ Foam::labelList Foam::surfaceFeatures::trimFeatures
 }
 
 
+void Foam::surfaceFeatures::excludeBox
+(
+    List<edgeStatus>& edgeStat,
+    const treeBoundBox& bb
+) const
+{
+    deleteBox(edgeStat, bb, true);
+}
+
+
+void Foam::surfaceFeatures::subsetBox
+(
+    List<edgeStatus>& edgeStat,
+    const treeBoundBox& bb
+) const
+{
+    deleteBox(edgeStat, bb, false);
+}
+
+
+void Foam::surfaceFeatures::deleteBox
+(
+    List<edgeStatus>& edgeStat,
+    const treeBoundBox& bb,
+    const bool removeInside
+) const
+{
+    const edgeList& surfEdges = surf_.edges();
+    const pointField& surfLocalPoints = surf_.localPoints();
+
+    forAll(edgeStat, edgei)
+    {
+        const point eMid = surfEdges[edgei].centre(surfLocalPoints);
+
+        if (removeInside ? bb.contains(eMid) : !bb.contains(eMid))
+        {
+            edgeStat[edgei] = surfaceFeatures::NONE;
+        }
+    }
+}
+
+
+void Foam::surfaceFeatures::subsetPlane
+(
+    List<edgeStatus>& edgeStat,
+    const plane& cutPlane
+) const
+{
+    const edgeList& surfEdges = surf_.edges();
+    const pointField& pts = surf_.points();
+    const labelList& meshPoints = surf_.meshPoints();
+
+    forAll(edgeStat, edgei)
+    {
+        const edge& e = surfEdges[edgei];
+
+        const point& p0 = pts[meshPoints[e.start()]];
+        const point& p1 = pts[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::surfaceFeatures::excludeOpen
+(
+    List<edgeStatus>& edgeStat
+) const
+{
+    forAll(edgeStat, edgei)
+    {
+        if (surf_.edgeFaces()[edgei].size() == 1)
+        {
+            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
+void Foam::surfaceFeatures::checkFlatRegionEdge
+(
+    List<edgeStatus>& edgeStat,
+    const scalar tol,
+    const scalar includedAngle
+) const
+{
+    forAll(edgeStat, edgei)
+    {
+        if (edgeStat[edgei] == surfaceFeatures::REGION)
+        {
+            const labelList& eFaces = surf_.edgeFaces()[edgei];
+
+            if (eFaces.size() > 2 && (eFaces.size() % 2) == 0)
+            {
+                edgeStat[edgei] = checkFlatRegionEdge
+                (
+                    tol,
+                    includedAngle,
+                    edgei
+                );
+            }
+        }
+    }
+}
+
+
 void Foam::surfaceFeatures::writeDict(Ostream& os) const
 {
-
     dictionary featInfoDict;
     featInfoDict.add("externalStart", externalStart_);
     featInfoDict.add("internalStart", internalStart_);
@@ -903,6 +1219,18 @@ void Foam::surfaceFeatures::writeObj(const fileName& prefix) const
 }
 
 
+void Foam::surfaceFeatures::writeStats(Ostream& os) const
+{
+    os  << "Feature set:" << nl
+        << "    points : " << this->featurePoints().size() << nl
+        << "    edges  : " << this->featureEdges().size() << nl
+        << "    of which" << nl
+        << "        region edges   : " << this->nRegionEdges() << nl
+        << "        external edges : " << this->nExternalEdges() << nl
+        << "        internal edges : " << this->nInternalEdges() << endl;
+}
+
+
 // Get nearest vertex on patch for every point of surf in pointSet.
 Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
 (
diff --git a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
index a9e12e267f..a4a3083f1a 100644
--- a/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
+++ b/src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2017 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -60,7 +60,9 @@ namespace Foam
 {
 
 // Forward declaration of classes
+class plane;
 class triSurface;
+class treeBoundBox;
 
 /*---------------------------------------------------------------------------*\
                            Class surfaceFeatures Declaration
@@ -73,10 +75,10 @@ public:
         //- Edge status
         enum edgeStatus
         {
-            NONE,
-            REGION,
-            EXTERNAL,
-            INTERNAL
+            NONE,       //!< Not a classified feature edge
+            REGION,     //
+            EXTERNAL,   //!< "Convex" edge
+            INTERNAL    //!< "Concave" edge
         };
 
 
@@ -170,6 +172,17 @@ private:
             labelList& featVisited
         );
 
+        //- Divide into multiple normal bins
+        //  - set REGION if != 2 normals
+        //  - set REGION if 2 normals that make feature angle
+        //  - otherwise set NONE and set normals,bins
+        edgeStatus checkFlatRegionEdge
+        (
+            const scalar tol,
+            const scalar includedAngle,
+            const label edgeI
+        ) const;
+
 public:
 
     ClassName("surfaceFeatures");
@@ -300,6 +313,50 @@ public:
                 const scalar includedAngle
             );
 
+            //- Mark edge status inside box as 'NONE'
+            void excludeBox
+            (
+                List<edgeStatus>& edgeStat,
+                const treeBoundBox& bb
+            ) const;
+
+            //- Mark edge status outside box as 'NONE'
+            void subsetBox
+            (
+                List<edgeStatus>& edgeStat,
+                const treeBoundBox& bb
+            ) const;
+
+            //- Mark edge status as 'NONE' for edges inside/outside box.
+            void deleteBox
+            (
+                List<edgeStatus>& edgeStat,
+                const treeBoundBox& bb,
+                const bool removeInside
+            ) const;
+
+            //- If edge does not intersect the plane, mark as 'NONE'
+            void subsetPlane
+            (
+                List<edgeStatus>& edgeStat,
+                const plane& cutPlane
+            ) const;
+
+            //- Mark edges with a single connected face as 'NONE'
+            void excludeOpen(List<edgeStatus>& edgeStat) const;
+
+            //- Divide into multiple normal bins
+            //  - set REGION if != 2 normals
+            //  - set REGION if 2 normals that make feature angle
+            //  - otherwise set NONE and set normals,bins
+            void checkFlatRegionEdge
+            (
+                List<edgeStatus>& edgeStat,
+                const scalar tol,
+                const scalar includedAngle
+            ) const;
+
+
             //- From member feature edges to status per edge.
             List<edgeStatus> toStatus() const;
 
@@ -414,6 +471,9 @@ public:
             //  feature points) for visualization
             void writeObj(const fileName& prefix) const;
 
+            //- Write some information
+            void writeStats(Ostream& os) const;
+
 
 
     // Member Operators
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceCloseness.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCloseness.C
new file mode 100644
index 0000000000..3925950b73
--- /dev/null
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCloseness.C
@@ -0,0 +1,360 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triSurfaceTools.H"
+
+#include "triSurface.H"
+#include "triSurfaceMesh.H"
+#include "triSurfaceFields.H"
+#include "MeshedSurface.H"
+#include "OFstream.H"
+#include "unitConversion.H"
+#include "meshTools.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+static 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;
+    }
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::Pair<Foam::tmp<Foam::scalarField>>
+Foam::triSurfaceTools::writeCloseness
+(
+    const Time& runTime,
+    const word& basename,
+    const triSurface& surf,
+    const scalar internalAngleTolerance,
+    const scalar externalAngleTolerance
+)
+{
+    Pair<tmp<scalarField>> tpair
+    (
+        tmp<scalarField>(new scalarField(surf.size(), GREAT)),
+        tmp<scalarField>(new scalarField(surf.size(), GREAT))
+    );
+
+    Info<< "Extracting internal and external closeness of surface." << endl;
+
+    triSurfaceMesh searchSurf
+    (
+        IOobject
+        (
+            basename + ".closeness",
+            runTime.constant(),
+            "triSurface",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf
+    );
+
+    // Prepare start and end points for intersection tests
+
+    const vectorField& normals = searchSurf.faceNormals();
+    const scalar span = searchSurf.bounds().mag();
+
+    const scalar externalToleranceCosAngle =
+        Foam::cos
+        (
+            degToRad(180 - externalAngleTolerance)
+        );
+
+    const scalar internalToleranceCosAngle =
+        Foam::cos
+        (
+            degToRad(180 - internalAngleTolerance)
+        );
+
+    Info<< "externalToleranceCosAngle: " << externalToleranceCosAngle << nl
+        << "internalToleranceCosAngle: " << internalToleranceCosAngle << endl;
+
+    // Info<< "span " << span << endl;
+
+    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 = tpair[0].ref();
+    scalarField& externalCloseness  = tpair[1].ref();
+
+    forAll(allHitInfo, fI)
+    {
+        const List<pointIndexHit>& hitInfo = allHitInfo[fI];
+
+        if (hitInfo.size() < 1)
+        {
+            drawHitProblem(fI, surf, start, faceCentres, end, hitInfo);
+
+            // FatalErrorInFunction
+            //     << "findLineAll did not hit its own face."
+            //     << exit(FatalError);
+        }
+        else if (hitInfo.size() == 1)
+        {
+            if (!hitInfo[0].hit())
+            {
+                // FatalErrorInFunction
+                //     << "findLineAll did not hit any face."
+                //     << exit(FatalError);
+            }
+            else if (hitInfo[0].index() != fI)
+            {
+                drawHitProblem
+                (
+                    fI,
+                    surf,
+                    start,
+                    faceCentres,
+                    end,
+                    hitInfo
+                );
+
+                // FatalErrorInFunction
+                //     << "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
+                );
+
+                // FatalErrorInFunction
+                //     << "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()
+                        );
+                }
+            }
+        }
+    }
+
+    // write as 'internalCloseness'
+    {
+        triSurfaceScalarField outputField
+        (
+            IOobject
+            (
+                basename + ".internalCloseness",
+                runTime.constant(),
+                "triSurface",
+                runTime,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            surf,
+            dimLength,
+            scalarField()
+        );
+
+        outputField.swap(internalCloseness);
+        outputField.write();
+        outputField.swap(internalCloseness);
+    }
+
+    // write as 'externalCloseness'
+    {
+        triSurfaceScalarField outputField
+        (
+            IOobject
+            (
+                basename + ".externalCloseness",
+                runTime.constant(),
+                "triSurface",
+                runTime,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            surf,
+            dimLength,
+            scalarField()
+        );
+
+        outputField.swap(externalCloseness);
+        outputField.write();
+        outputField.swap(externalCloseness);
+    }
+
+    return tpair;
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C
new file mode 100644
index 0000000000..2518b09516
--- /dev/null
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceCurvature.C
@@ -0,0 +1,386 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triSurfaceTools.H"
+
+#include "triSurface.H"
+#include "MeshedSurfaces.H"
+#include "triSurfaceFields.H"
+#include "OFstream.H"
+#include "plane.H"
+#include "tensor2D.H"
+#include "symmTensor2D.H"
+#include "scalarMatrices.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::triSurfaceTools::vertexNormalWeight
+(
+    const triFace& f,
+    const label pI,
+    const vector& fN,
+    const UList<point>& 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::tmp<Foam::vectorField>
+Foam::triSurfaceTools::vertexNormals(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;
+
+    tmp<vectorField> tfld(new vectorField(surf.nPoints(), Zero));
+    vectorField& pointNormals = tfld.ref();
+
+    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);
+
+            const scalar weight = vertexNormalWeight
+            (
+                f,
+                meshPoints[pI],
+                fN,
+                points
+            );
+
+            pointNormals[pI] += weight*fN;
+        }
+
+        pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL;
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::triadField>
+Foam::triSurfaceTools::vertexTriads
+(
+    const triSurface& surf,
+    const vectorField& pointNormals
+)
+{
+    const pointField& points = surf.points();
+    const Map<label>& meshPointMap = surf.meshPointMap();
+
+    tmp<triadField> tfld(new triadField(points.size()));
+    triadField& pointTriads = tfld.ref();
+
+    forAll(points, pI)
+    {
+        const point& pt = points[pI];
+        const vector& normal = pointNormals[meshPointMap[pI]];
+
+        if (mag(normal) < SMALL)
+        {
+            pointTriads[meshPointMap[pI]] = triad::unset;
+            continue;
+        }
+
+        plane p(pt, normal);
+
+        // Pick arbitrary point in plane
+        vector dir1 = pt - p.somePointInPlane(1e-3);
+        dir1 /= mag(dir1);
+
+        vector dir2 = dir1 ^ normal;
+        dir2 /= mag(dir2);
+
+        pointTriads[meshPointMap[pI]] = triad(dir1, dir2, normal);
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::triSurfaceTools::curvatures
+(
+    const triSurface& surf,
+    const vectorField& pointNormals,
+    const triadField& pointTriads
+)
+{
+    Info<< "Calculating face curvature" << endl;
+
+    const pointField& points = surf.points();
+    const labelList& meshPoints = surf.meshPoints();
+    const Map<label>& meshPointMap = surf.meshPointMap();
+
+    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 = pointTriads[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 = triSurfaceTools::vertexNormalWeight
+            (
+                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;
+        }
+    }
+
+    tmp<scalarField> tfld(new scalarField(points.size(), Zero));
+    scalarField& curvatureAtPoints = tfld.ref();
+
+    forAll(curvatureAtPoints, 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];
+
+        curvatureAtPoints[meshPoints[pI]] = curvature;
+    }
+
+    return tfld;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::triSurfaceTools::curvatures
+(
+    const triSurface& surf
+)
+{
+    tmp<vectorField> norms = triSurfaceTools::vertexNormals(surf);
+    tmp<triadField> triads = triSurfaceTools::vertexTriads(surf, norms());
+
+    tmp<scalarField> curv = curvatures(surf, norms(), triads());
+    norms.clear();
+    triads.clear();
+
+    return curv;
+}
+
+
+Foam::tmp<Foam::scalarField>
+Foam::triSurfaceTools::writeCurvature
+(
+    const Time& runTime,
+    const word& basename,
+    const triSurface& surf
+)
+{
+    Info<< "Extracting curvature of surface at the points." << endl;
+
+    tmp<scalarField> tfld = triSurfaceTools::curvatures(surf);
+    scalarField& curv = tfld.ref();
+
+    triSurfacePointScalarField outputField
+    (
+        IOobject
+        (
+            basename + ".curvature",
+            runTime.constant(),
+            "triSurface",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        scalarField()
+    );
+
+    outputField.swap(curv);
+    outputField.write();
+    outputField.swap(curv);
+
+    return tfld;
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index aec02f8e06..d1da45a001 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -33,7 +33,6 @@ License
 #include "plane.H"
 #include "geompack.H"
 
-
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const Foam::label Foam::triSurfaceTools::ANYEDGE = -1;
@@ -1376,54 +1375,6 @@ Foam::labelList Foam::triSurfaceTools::getVertexVertices
 }
 
 
-//// Order vertices consistent with face
-//void Foam::triSurfaceTools::orderVertices
-//(
-//    const labelledTri& f,
-//    const label v1,
-//    const label v2,
-//    label& vA,
-//    label& vB
-//)
-//{
-//    // Order v1, v2 in anticlockwise order.
-//    bool reverse = false;
-//
-//    if (f[0] == v1)
-//    {
-//        if (f[1] != v2)
-//        {
-//            reverse = true;
-//        }
-//    }
-//    else if (f[1] == v1)
-//    {
-//        if (f[2] != v2)
-//        {
-//            reverse = true;
-//        }
-//    }
-//    else
-//    {
-//        if (f[0] != v2)
-//        {
-//            reverse = true;
-//        }
-//    }
-//
-//    if (reverse)
-//    {
-//        vA = v2;
-//        vB = v1;
-//    }
-//    else
-//    {
-//        vA = v1;
-//        vB = v2;
-//    }
-//}
-
-
 // Get the other face using edgeI
 Foam::label Foam::triSurfaceTools::otherFace
 (
@@ -2633,7 +2584,6 @@ void Foam::triSurfaceTools::calcInterpolationWeights
     {
         const point& samplePt = samplePts[i];
 
-
         FixedList<label, 3>& verts = allVerts[i];
         FixedList<scalar, 3>& weights = allWeights[i];
 
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
index d90c862dc9..28960823db 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
@@ -27,8 +27,21 @@ Class
 Description
     A collection of tools for triSurface.
 
+Note
+    The curvature calculation is an implementation of the algorithm from:
+    \verbatim
+        "Estimating Curvatures and their Derivatives on Triangle Meshes"
+        by S. Rusinkiewicz
+        3DPVT'04 Proceedings of the 3D Data Processing,
+        Visualization, and Transmission, 2nd International Symposium
+        Pages 486-493
+        http://gfx.cs.princeton.edu/pubs/_2004_ECA/curvpaper.pdf
+    \endverbatim
+
 SourceFiles
     triSurfaceTools.C
+    triSurfaceCloseness.C
+    triSurfaceCurvature.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -37,9 +50,12 @@ SourceFiles
 
 #include "boolList.H"
 #include "pointField.H"
+#include "vectorField.H"
+#include "triadFieldFwd.H"
 #include "DynamicList.H"
 #include "HashSet.H"
 #include "FixedList.H"
+#include "Pair.H"
 #include "vector2D.H"
 #include "triPointRef.H"
 #include "surfaceLocation.H"
@@ -56,6 +72,7 @@ class polyBoundaryMesh;
 class plane;
 class triSurface;
 class face;
+class Time;
 template<class Face> class MeshedSurface;
 
 
@@ -280,16 +297,6 @@ public:
             const edge& e
         );
 
-        ////- Order vertices consistent with face
-        //static void orderVertices
-        //(
-        //    const labelledTri& f,
-        //    const label v1,
-        //    const label v2,
-        //    label& vA,
-        //    label& vB
-        //);
-
         //- Get face connected to edge not facei
         static label otherFace
         (
@@ -490,12 +497,16 @@ public:
 
     // Triangulation and interpolation
 
+        //- Do unconstrained Delaunay of points. Returns triSurface with 3D
+        //  points with z=0. All triangles in region 0.
+        static triSurface delaunay2D(const List<vector2D>&);
+
         //- Calculate linear interpolation weights for point (guaranteed to be
         //  inside triangle)
         static void calcInterpolationWeights
         (
-            const triPointRef&,
-            const point&,
+            const triPointRef& tri,
+            const point& p,
             FixedList<scalar, 3>& weights
         );
 
@@ -510,13 +521,86 @@ public:
         (
             const triSurface& s,
             const pointField& samplePts,
-            List<FixedList<label, 3>>& verts,
-            List<FixedList<scalar, 3>>& weights
+            List<FixedList<label, 3>>& allVerts,
+            List<FixedList<scalar, 3>>& allWeights
         );
 
-        //- Do unconstrained Delaunay of points. Returns triSurface with 3D
-        //  points with z=0. All triangles in region 0.
-        static triSurface delaunay2D(const List<vector2D>&);
+
+    // Curvature
+
+        //- Weighting for normals of faces attached to vertex
+        static scalar vertexNormalWeight
+        (
+            const triFace& f,
+            const label pI,
+            const vector& fN,
+            const UList<point>& points
+        );
+
+        //- Weighted average of normals of attached faces
+        static tmp<vectorField> vertexNormals(const triSurface& surf);
+
+        //- Local coordinate-system for each point normal
+        static tmp<triadField> vertexTriads
+        (
+            const triSurface& surf,
+            const vectorField& pointNormals
+        );
+
+        //- Surface curvatures at the vertex points
+        static tmp<scalarField> curvatures
+        (
+            const triSurface& surf,
+            const vectorField& pointNormals,
+            const triadField& pointTriads
+        );
+
+        //- Surface curvatures at the vertex points
+        static tmp<scalarField> curvatures
+        (
+            const triSurface& surf
+        );
+
+        //- Write surface curvature at the vertex points and return the field
+        static tmp<scalarField> writeCurvature
+        (
+            const Time& runTime,
+            const word& basename,
+            const triSurface& surf
+        );
+
+
+    // Closeness
+
+        //- Check and write internal/external closeness fields
+        static Pair<tmp<scalarField>> writeCloseness
+        (
+            const Time& runTime,
+            const word& basename,
+            const triSurface& surf,
+            const scalar internalAngleTolerance = 45,
+            const scalar externalAngleTolerance = 10
+        );
+
+
+    // Feature Proximity
+
+        //- Calculate feature proximity
+        scalarField featureProximity
+        (
+            const triSurface& surf,
+            const scalar searchDistance
+        );
+
+        //- Check and write internal/external closeness fields
+        static void writeFeatureProximity
+        (
+            const Time& runTime,
+            const word& basename,
+            const triSurface& surf,
+            const bool writeVTK,
+            const scalar searchDistance
+        );
 
 
     // Surface checking functionality