ENH: refactor code from surfaceFeatureExtract to triSurfaceTools (issue #450)

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
 

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;
+}
+
+// ************************************************************************* //

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;
+}
+
+
+// ************************************************************************* //

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
+
+// ************************************************************************* //

src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.H

@@ -94,23 +94,22 @@ public:
 
     enum pointStatus
     {
-        CONVEX,     // Fully convex point (w.r.t normals)
+        CONVEX,     //!< Fully convex point (w.r.t normals)
-        CONCAVE,    // Fully concave point
+        CONCAVE,    //!< Fully concave point
-        MIXED,      // A point surrounded by both convex and concave edges
+        MIXED,      //!< A point surrounded by both convex and concave edges
-        NONFEATURE  // Not a feature point
+        NONFEATURE  //!< Not a feature point
     };
 
     static const Foam::NamedEnum<pointStatus, 4> pointStatusNames_;
 
     enum edgeStatus
     {
-        EXTERNAL,   // "Convex" edge
+        EXTERNAL,   //!< "Convex" edge
-        INTERNAL,   // "Concave" edge
+        INTERNAL,   //!< "Concave" edge
-        FLAT,       // Neither concave or convex, on a flat surface
+        FLAT,       //!< Neither concave or convex, on a flat surface
-        OPEN,       // i.e. only connected to one face
+        OPEN,       //!< Only connected to a single face
-        MULTIPLE,   // Multiply connected (connected to more than two faces)
+        MULTIPLE,   //!< Multiply connected (connected to more than two faces)
-        NONE        // Not a classified feature edge (consistency with
+        NONE        //!< Unclassified (consistency with surfaceFeatures)
-                    // surfaceFeatures)
     };
 
     static const Foam::NamedEnum<edgeStatus, 6> edgeStatusNames_;
@@ -118,10 +117,10 @@ public:
     //- Normals point to the outside
     enum sideVolumeType
     {
-        INSIDE  = 0,  // mesh inside
+        INSIDE  = 0,  //!< mesh inside
-        OUTSIDE = 1,  // mesh outside
+        OUTSIDE = 1,  //!< mesh outside
-        BOTH    = 2,  // e.g. a baffle
+        BOTH    = 2,  //!< e.g. a baffle
-        NEITHER = 3   // not sure when this may be used
+        NEITHER = 3   //!< not sure when this may be used
     };
 
     static const Foam::NamedEnum<sideVolumeType, 4> sideVolumeTypeNames_;

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
 (

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,
+            NONE,       //!< Not a classified feature edge
-            REGION,
+            REGION,     //
-            EXTERNAL,
+            EXTERNAL,   //!< "Convex" edge
-            INTERNAL
+            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

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;
+}
+
+
+// ************************************************************************* //

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;
+}
+
+
+// ************************************************************************* //

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];
 

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 triPointRef& tri,
-            const point&,
+            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<label, 3>>& allVerts,
-            List<FixedList<scalar, 3>>& weights
+            List<FixedList<scalar, 3>>& allWeights
         );
 
-        //- Do unconstrained Delaunay of points. Returns triSurface with 3D
+
-        //  points with z=0. All triangles in region 0.
+    // Curvature
-        static triSurface delaunay2D(const List<vector2D>&);
+
+        //- 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