From 21d5267b6fb67d4fdcb4461a0c3998ce3cf5dec6 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Fri, 6 Apr 2018 20:15:28 +0100
Subject: [PATCH] surfaceFeatureExtract: Refactor core functionality into core
 classes and libraries

---
 .../surfaceFeatureExtract.C                   |  30 +-
 .../surfaceFeatureExtract.H                   |  17 -
 .../surfaceFeatureExtractUtilities.C          | 335 ----------
 .../meshes/primitiveShapes/plane/plane.C      |  67 ++
 .../meshes/primitiveShapes/plane/plane.H      |   5 +-
 .../meshes/treeBoundBox/treeBoundBox.C        |  26 +
 .../meshes/treeBoundBox/treeBoundBox.H        |   6 +
 src/triSurface/triSurface/triSurface.C        | 590 ++++++++++++------
 src/triSurface/triSurface/triSurface.H        |  33 +-
 9 files changed, 533 insertions(+), 576 deletions(-)

diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
index 60f682529..fdb3c9ec3 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C
@@ -231,18 +231,18 @@ int main(int argc, char *argv[])
             {
                 treeBoundBox bb(subsetDict.lookup("insideBox")());
 
-                Info<< "Removing all edges outside bb " << bb << endl;
-                dumpBox(bb, "subsetBox.obj");
-
+                Info<< "Removing all edges outside bb " << bb
+                    << " see subsetBox.obj" << endl;
+                bb.writeOBJ("subsetBox.obj");
                 deleteBox(surf, bb, false, edgeStat);
             }
             else if (subsetDict.found("outsideBox"))
             {
                 treeBoundBox bb(subsetDict.lookup("outsideBox")());
 
-                Info<< "Removing all edges inside bb " << bb << endl;
-                dumpBox(bb, "deleteBox.obj");
-
+                Info<< "Removing all edges inside bb " << bb
+                    << " see deleteBox.obj" << endl;
+                bb.writeOBJ("deleteBox.obj");
                 deleteBox(surf, bb, true, edgeStat);
             }
 
@@ -432,16 +432,18 @@ int main(int argc, char *argv[])
             Info<< nl << "Extracting curvature of surface at the points."
                 << endl;
 
-            const vectorField pointNormals(surf.pointNormals2());
-            triadField pointCoordSys = calcVertexCoordSys(surf, pointNormals);
-
-            triSurfacePointScalarField k = calcCurvature
+            triSurfacePointScalarField k
             (
-                sFeatFileName,
-                runTime,
+                IOobject
+                (
+                    sFeatFileName + ".curvature",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime
+                ),
                 surf,
-                pointNormals,
-                pointCoordSys
+                dimLength,
+                surf.curvature()
             );
 
             k.write();
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H
index 40333d3c1..7c326975e 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H
@@ -45,21 +45,6 @@ namespace Foam
 
     point randomPointInPlane(const plane& p);
 
-    triadField calcVertexCoordSys
-    (
-        const triSurface& surf,
-        const vectorField& pointNormals
-    );
-
-    triSurfacePointScalarField calcCurvature
-    (
-        const word& name,
-        const Time& runTime,
-        const triSurface& surf,
-        const vectorField& pointNormals,
-        const triadField& pointCoordSys
-    );
-
     bool edgesConnected(const edge& e1, const edge& e2);
 
     scalar calcProximityOfFeaturePoints
@@ -75,8 +60,6 @@ namespace Foam
         const scalar defaultCellSize
     );
 
-    void dumpBox(const treeBoundBox& bb, const fileName& fName);
-
     //- Deletes all edges inside/outside bounding box from set.
     void deleteBox
     (
diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C
index 2321fb8b6..42a5e083c 100644
--- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C
+++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C
@@ -28,11 +28,6 @@ Description
     Extracts and writes surface features to file. All but the basic feature
     extraction is WIP.
 
-    Curvature calculation is an implementation of the algorithm from:
-
-      "Estimating Curvatures and their Derivatives on Triangle Meshes"
-      by S. Rusinkiewicz
-
 \*---------------------------------------------------------------------------*/
 
 #include "surfaceFeatureExtract.H"
@@ -56,294 +51,6 @@ const Foam::scalar Foam::externalToleranceCosAngle
 );
 
 
-Foam::point Foam::randomPointInPlane(const plane& p)
-{
-    // Perturb base point
-    const point& refPt = p.refPoint();
-
-    // ax + by + cz + d = 0
-    const FixedList<scalar, 4>& planeCoeffs = p.planeCoeffs();
-
-    const scalar perturbX = refPt.x() + 1e-3;
-    const scalar perturbY = refPt.y() + 1e-3;
-    const scalar perturbZ = refPt.z() + 1e-3;
-
-    if (mag(planeCoeffs[2]) < small)
-    {
-        if (mag(planeCoeffs[1]) < small)
-        {
-            const scalar x =
-                -1.0
-                *(
-                     planeCoeffs[3]
-                   + planeCoeffs[1]*perturbY
-                   + planeCoeffs[2]*perturbZ
-                 )/planeCoeffs[0];
-
-            return point
-            (
-                x,
-                perturbY,
-                perturbZ
-            );
-        }
-
-        const scalar y =
-            -1.0
-            *(
-                 planeCoeffs[3]
-               + planeCoeffs[0]*perturbX
-               + planeCoeffs[2]*perturbZ
-             )/planeCoeffs[1];
-
-        return point
-        (
-            perturbX,
-            y,
-            perturbZ
-        );
-    }
-    else
-    {
-        const scalar z =
-            -1.0
-            *(
-                 planeCoeffs[3]
-               + planeCoeffs[0]*perturbX
-               + planeCoeffs[1]*perturbY
-             )/planeCoeffs[2];
-
-        return point
-        (
-            perturbX,
-            perturbY,
-            z
-        );
-    }
-}
-
-
-Foam::triadField Foam::calcVertexCoordSys
-(
-    const triSurface& surf,
-    const vectorField& pointNormals
-)
-{
-    const pointField& points = surf.points();
-    const Map<label>& meshPointMap = surf.meshPointMap();
-
-    triadField pointCoordSys(points.size());
-
-    forAll(points, pI)
-    {
-        const point& pt = points[pI];
-        const vector& normal = pointNormals[meshPointMap[pI]];
-
-        if (mag(normal) < small)
-        {
-            pointCoordSys[meshPointMap[pI]] = triad::unset;
-            continue;
-        }
-
-        plane p(pt, normal);
-
-        // Pick random point in plane
-        vector dir1 = pt - randomPointInPlane(p);
-        dir1 /= mag(dir1);
-
-        vector dir2 = dir1 ^ normal;
-        dir2 /= mag(dir2);
-
-        pointCoordSys[meshPointMap[pI]] = triad(dir1, dir2, normal);
-    }
-
-    return pointCoordSys;
-}
-
-
-Foam::triSurfacePointScalarField Foam::calcCurvature
-(
-    const word& name,
-    const Time& runTime,
-    const triSurface& surf,
-    const vectorField& pointNormals,
-    const triadField& pointCoordSys
-)
-{
-    Info<< "Calculating face curvature" << endl;
-
-    const pointField& points = surf.points();
-    const labelList& meshPoints = surf.meshPoints();
-    const Map<label>& meshPointMap = surf.meshPointMap();
-
-    triSurfacePointScalarField curvaturePointField
-    (
-        IOobject
-        (
-            name + ".curvature",
-            runTime.constant(),
-            "triSurface",
-            runTime,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        surf,
-        dimLength,
-        scalarField(points.size(), 0.0)
-    );
-
-    List<symmTensor2D> pointFundamentalTensors
-    (
-        points.size(),
-        symmTensor2D::zero
-    );
-
-    scalarList accumulatedWeights(points.size(), 0.0);
-
-    forAll(surf, fi)
-    {
-        const triFace& f = surf[fi];
-        const edgeList fEdges = f.edges();
-
-        // Calculate the edge vectors and the normal differences
-        vectorField edgeVectors(f.size(), Zero);
-        vectorField normalDifferences(f.size(), Zero);
-
-        forAll(fEdges, feI)
-        {
-            const edge& e = fEdges[feI];
-
-            edgeVectors[feI] = e.vec(points);
-            normalDifferences[feI] =
-               pointNormals[meshPointMap[e[0]]]
-             - pointNormals[meshPointMap[e[1]]];
-        }
-
-        // Set up a local coordinate system for the face
-        const vector& e0 = edgeVectors[0];
-        const vector eN = f.area(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
-            const scalar weight = surf.pointNormalWeight
-            (
-                f,
-                meshPoints[patchPointIndex],
-                f.area(points),
-                points
-            );
-
-            // Sum contribution of face to this point
-            pointFundamentalTensors[patchPointIndex] +=
-                weight*projectedFundamentalTensor;
-
-            accumulatedWeights[patchPointIndex] += weight;
-        }
-    }
-
-    forAll(curvaturePointField, pI)
-    {
-        pointFundamentalTensors[pI] /= (accumulatedWeights[pI] + small);
-
-        vector2D principalCurvatures = eigenValues(pointFundamentalTensors[pI]);
-
-        //scalar curvature =
-        //    (principalCurvatures[0] + principalCurvatures[1])/2;
-        scalar curvature = max
-        (
-            mag(principalCurvatures[0]),
-            mag(principalCurvatures[1])
-        );
-        //scalar curvature = principalCurvatures[0]*principalCurvatures[1];
-
-        curvaturePointField[meshPoints[pI]] = curvature;
-    }
-
-    return curvaturePointField;
-}
-
-
 bool Foam::edgesConnected(const edge& e1, const edge& e2)
 {
     if
@@ -455,30 +162,6 @@ Foam::scalar Foam::calcProximityOfFeatureEdges
 }
 
 
-void Foam::dumpBox(const treeBoundBox& bb, const fileName& fName)
-{
-    OFstream str(fName);
-
-    Info<< "Dumping bounding box " << bb << " as lines to obj file "
-        << str.name() << endl;
-
-
-    pointField boxPoints(bb.points());
-
-    forAll(boxPoints, i)
-    {
-        meshTools::writeOBJ(str, boxPoints[i]);
-    }
-
-    forAll(treeBoundBox::edges, i)
-    {
-        const edge& e = treeBoundBox::edges[i];
-
-        str<< "l " << e[0]+1 <<  ' ' << e[1]+1 << nl;
-    }
-}
-
-
 void Foam::deleteBox
 (
     const triSurface& surf,
@@ -809,24 +492,6 @@ Foam::surfaceFeatures::edgeStatus Foam::checkFlatRegionEdge
 }
 
 
-void Foam::extractCloseness
-(
-    const fileName &sFeatFileName,
-    const Time& runTime,
-    const triSurface &surf,
-    const bool writeVTK
-);
-
-
-void Foam::extractPointCloseness
-(
-    const fileName &sFeatFileName,
-    const Time& runTime,
-    const triSurface &surf,
-    const bool writeVTK
-);
-
-
 void Foam::writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
 {
     os  << "    points : " << fem.points().size() << nl
diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
index 81408d164..f89ded9cc 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C
@@ -289,6 +289,73 @@ Foam::FixedList<Foam::scalar, 4> Foam::plane::planeCoeffs() const
 }
 
 
+Foam::point Foam::plane::aPoint() const
+{
+    // Perturb base point
+    const point& refPt = refPoint();
+
+    // ax + by + cz + d = 0
+    FixedList<scalar, 4> planeCoeffs = this->planeCoeffs();
+
+    const scalar perturbX = refPt.x() + 1e-3;
+    const scalar perturbY = refPt.y() + 1e-3;
+    const scalar perturbZ = refPt.z() + 1e-3;
+
+    if (mag(planeCoeffs[2]) < small)
+    {
+        if (mag(planeCoeffs[1]) < small)
+        {
+            const scalar x =
+                -1.0
+                *(
+                     planeCoeffs[3]
+                   + planeCoeffs[1]*perturbY
+                   + planeCoeffs[2]*perturbZ
+                 )/planeCoeffs[0];
+
+            return point
+            (
+                x,
+                perturbY,
+                perturbZ
+            );
+        }
+
+        const scalar y =
+            -1.0
+            *(
+                 planeCoeffs[3]
+               + planeCoeffs[0]*perturbX
+               + planeCoeffs[2]*perturbZ
+             )/planeCoeffs[1];
+
+        return point
+        (
+            perturbX,
+            y,
+            perturbZ
+        );
+    }
+    else
+    {
+        const scalar z =
+            -1.0
+            *(
+                 planeCoeffs[3]
+               + planeCoeffs[0]*perturbX
+               + planeCoeffs[1]*perturbY
+             )/planeCoeffs[2];
+
+        return point
+        (
+            perturbX,
+            perturbY,
+            z
+        );
+    }
+}
+
+
 Foam::point Foam::plane::nearestPoint(const point& p) const
 {
     return p - normal_*((p - point_) & normal_);
diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
index 5265ef063..a934d1ae0 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -162,6 +162,9 @@ public:
         //  plane equation: ax + by + cz + d = 0
         FixedList<scalar, 4> planeCoeffs() const;
 
+        //- Return a point on the plane
+        point aPoint() const;
+
         //- Return nearest point in the plane for the given point
         point nearestPoint(const point& p) const;
 
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
index a38a46040..0eb41a1ef 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
+++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C
@@ -25,6 +25,7 @@ License
 
 #include "treeBoundBox.H"
 #include "ListOps.H"
+#include "OFstream.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -615,6 +616,31 @@ Foam::label Foam::treeBoundBox::distanceCmp
 }
 
 
+void Foam::treeBoundBox::writeOBJ(const fileName& fName) const
+{
+    OFstream str(fName);
+
+    Info<< "Dumping bounding box " << *this << " as lines to obj file "
+        << str.name() << endl;
+
+    pointField bbPoints(points());
+
+    forAll(bbPoints, i)
+    {
+        const point& pt = bbPoints[i];
+
+        str<< "v " << pt.x() << ' ' << pt.y() << ' ' << pt.z() << endl;
+    }
+
+    forAll(treeBoundBox::edges, i)
+    {
+        const edge& e = treeBoundBox::edges[i];
+
+        str<< "l " << e[0]+1 <<  ' ' << e[1]+1 << nl;
+    }
+}
+
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 bool Foam::operator==(const treeBoundBox& a, const treeBoundBox& b)
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
index 955b01e42..8c35577fd 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
+++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
@@ -349,11 +349,17 @@ public:
             //  and guarantees in any direction s*mag(span) minimum width
             inline treeBoundBox extend(Random&, const scalar s) const;
 
+    // Write
+
+        void writeOBJ(const fileName& fName) const;
+
+
     // Friend Operators
 
         friend bool operator==(const treeBoundBox&, const treeBoundBox&);
         friend bool operator!=(const treeBoundBox&, const treeBoundBox&);
 
+
     // IOstream operator
 
         friend Istream& operator>>(Istream& is, treeBoundBox&);
diff --git a/src/triSurface/triSurface/triSurface.C b/src/triSurface/triSurface/triSurface.C
index da72f5511..90aa10d92 100644
--- a/src/triSurface/triSurface/triSurface.C
+++ b/src/triSurface/triSurface/triSurface.C
@@ -31,6 +31,10 @@ License
 #include "boundBox.H"
 #include "SortableList.H"
 #include "PackedBoolList.H"
+#include "plane.H"
+#include "tensor2D.H"
+#include "symmTensor2D.H"
+#include "transform.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -181,164 +185,105 @@ Foam::string Foam::triSurface::getLineNoComment(IFstream& is)
 }
 
 
-// Remove non-triangles, double triangles.
-void Foam::triSurface::checkTriangles(const bool verbose)
+Foam::scalar Foam::triSurface::pointNormalWeight
+(
+    const triFace& f,
+    const label pi,
+    const vector& fa,
+    const pointField& points
+) const
 {
-    // Simple check on indices ok.
-    const label maxPointi = points().size() - 1;
+    const label index = findIndex(f, pi);
 
-    forAll(*this, facei)
+    if (index == -1)
     {
-        const triSurface::FaceType& f = (*this)[facei];
-
-        forAll(f, fp)
-        {
-            if (f[fp] < 0 || f[fp] > maxPointi)
-            {
-                FatalErrorInFunction
-                    << "triangle " << f
-                    << " uses point indices outside point range 0.."
-                    << maxPointi
-                    << exit(FatalError);
-            }
-        }
+        FatalErrorInFunction
+            << "Point not in face" << abort(FatalError);
     }
 
-    // Two phase process
-    //   1. mark invalid faces
-    //   2. pack
-    // Done to keep numbering constant in phase 1
+    const vector e1 = points[f[index]] - points[f[f.fcIndex(index)]];
+    const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
 
-    // List of valid triangles
-    boolList valid(size(), true);
-    bool hasInvalid = false;
-
-    forAll(*this, facei)
-    {
-        const labelledTri& f = (*this)[facei];
-
-        if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
-        {
-            // 'degenerate' triangle check
-            valid[facei] = false;
-            hasInvalid = true;
-
-            if (verbose)
-            {
-                WarningInFunction
-                    << "triangle " << facei
-                    << " does not have three unique vertices:\n";
-                printTriangle(Warning, "    ", f, points());
-            }
-        }
-        else
-        {
-            // duplicate triangle check
-            const labelList& fEdges = faceEdges()[facei];
-
-            // Check if faceNeighbours use same points as this face.
-            // Note: discards normal information - sides of baffle are merged.
-
-            forAll(fEdges, fp)
-            {
-                const labelList& eFaces = edgeFaces()[fEdges[fp]];
-
-                forAll(eFaces, i)
-                {
-                    label neighbour = eFaces[i];
-
-                    if (neighbour > facei)
-                    {
-                        // lower numbered faces already checked
-                        const labelledTri& n = (*this)[neighbour];
-
-                        if
-                        (
-                            ((f[0] == n[0]) || (f[0] == n[1]) || (f[0] == n[2]))
-                         && ((f[1] == n[0]) || (f[1] == n[1]) || (f[1] == n[2]))
-                         && ((f[2] == n[0]) || (f[2] == n[1]) || (f[2] == n[2]))
-                        )
-                        {
-                            valid[facei] = false;
-                            hasInvalid = true;
-
-                            if (verbose)
-                            {
-                                WarningInFunction
-                                    << "triangles share the same vertices:\n"
-                                    << "    face 1 :" << facei << endl;
-                                printTriangle(Warning, "    ", f, points());
-
-                                Warning
-                                    << endl
-                                    << "    face 2 :"
-                                    << neighbour << endl;
-                                printTriangle(Warning, "    ", n, points());
-                            }
-
-                            break;
-                        }
-                    }
-                }
-            }
-        }
-    }
-
-    if (hasInvalid)
-    {
-        // Pack
-        label newFacei = 0;
-        forAll(*this, facei)
-        {
-            if (valid[facei])
-            {
-                const labelledTri& f = (*this)[facei];
-                (*this)[newFacei++] = f;
-            }
-        }
-
-        if (verbose)
-        {
-            WarningInFunction
-                << "Removing " << size() - newFacei
-                << " illegal faces." << endl;
-        }
-        (*this).setSize(newFacei);
-
-        // Topology can change because of renumbering
-        clearOut();
-    }
+    return mag(fa)/(magSqr(e1)*magSqr(e2) + vSmall);
 }
 
 
-// Check/fix edges with more than two triangles
-void Foam::triSurface::checkEdges(const bool verbose)
+Foam::tmp<Foam::vectorField> Foam::triSurface::weightedPointNormals() const
 {
-    const labelListList& eFaces = edgeFaces();
+    // Weighted average of normals of faces attached to the vertex
+    // Weight = fA / (mag(e0)^2 * mag(e1)^2);
 
-    forAll(eFaces, edgeI)
+    tmp<vectorField> tpointNormals
+    (
+        new vectorField(nPoints(), Zero)
+    );
+    vectorField& pointNormals = tpointNormals.ref();
+
+    const pointField& points = this->points();
+    const labelListList& pointFaces = this->pointFaces();
+    const labelList& meshPoints = this->meshPoints();
+
+    forAll(pointFaces, pi)
     {
-        const labelList& myFaces = eFaces[edgeI];
+        const labelList& pFaces = pointFaces[pi];
 
-        if (myFaces.empty())
+        forAll(pFaces, fi)
         {
-            FatalErrorInFunction
-                << "Edge " << edgeI << " with vertices " << edges()[edgeI]
-                << " has no edgeFaces"
-                << exit(FatalError);
-        }
-        else if (myFaces.size() > 2 && verbose)
-        {
-            WarningInFunction
-                << "Edge " << edgeI << " with vertices " << edges()[edgeI]
-                << " has more than 2 faces connected to it : " << myFaces
-                << endl;
+            const label facei = pFaces[fi];
+            const triFace& f = operator[](facei);
+
+            const vector fa = f.area(points);
+            const scalar weight =
+                pointNormalWeight(f, meshPoints[pi], fa, points);
+
+            pointNormals[pi] += weight*fa;
         }
+
+        pointNormals[pi] /= mag(pointNormals[pi]) + vSmall;
     }
+
+    return tpointNormals;
+}
+
+
+Foam::tmp<Foam::triadField> Foam::triSurface::pointCoordSys
+(
+    const vectorField& pointNormals
+) const
+{
+    tmp<triadField> tpointCoordSys(new triadField(nPoints()));
+    triadField& pointCoordSys = tpointCoordSys.ref();
+
+    const pointField& points = this->points();
+    const labelList& meshPoints = this->meshPoints();
+
+    forAll(meshPoints, pi)
+    {
+        const point& pt = points[meshPoints[pi]];
+        const vector& normal = pointNormals[pi];
+
+        if (mag(normal) < small)
+        {
+            pointCoordSys[pi] = triad::unset;
+            continue;
+        }
+
+        plane p(pt, normal);
+
+        // Pick a point in plane
+        vector dir1 = pt - p.aPoint();
+        dir1 /= mag(dir1);
+
+        vector dir2 = dir1 ^ normal;
+        dir2 /= mag(dir2);
+
+        pointCoordSys[pi] = triad(dir1, dir2, normal);
+    }
+
+    return pointCoordSys;
 }
 
 
-// Read triangles, points from Istream
 bool Foam::triSurface::read(Istream& is)
 {
     is  >> patches_ >> storedPoints() >> storedFaces();
@@ -347,7 +292,6 @@ bool Foam::triSurface::read(Istream& is)
 }
 
 
-// Read from file in given format
 bool Foam::triSurface::read
 (
     const fileName& name,
@@ -421,7 +365,6 @@ bool Foam::triSurface::read
 }
 
 
-// Write to file in given format
 void Foam::triSurface::write
 (
     const fileName& name,
@@ -484,8 +427,6 @@ void Foam::triSurface::write
 }
 
 
-// Returns patch info. Sets faceMap to the indexing according to patch
-// numbers. Patch numbers start at 0.
 Foam::surfacePatchList Foam::triSurface::calcPatches(labelList& faceMap) const
 {
     // Sort according to region numbers of labelledTri
@@ -812,7 +753,161 @@ void Foam::triSurface::scalePoints(const scalar scaleFactor)
 }
 
 
-// Remove non-triangles, double triangles.
+void Foam::triSurface::checkTriangles(const bool verbose)
+{
+    // Simple check on indices ok.
+    const label maxPointi = points().size() - 1;
+
+    forAll(*this, facei)
+    {
+        const triSurface::FaceType& f = (*this)[facei];
+
+        forAll(f, fp)
+        {
+            if (f[fp] < 0 || f[fp] > maxPointi)
+            {
+                FatalErrorInFunction
+                    << "triangle " << f
+                    << " uses point indices outside point range 0.."
+                    << maxPointi
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    // Two phase process
+    //   1. mark invalid faces
+    //   2. pack
+    // Done to keep numbering constant in phase 1
+
+    // List of valid triangles
+    boolList valid(size(), true);
+    bool hasInvalid = false;
+
+    forAll(*this, facei)
+    {
+        const labelledTri& f = (*this)[facei];
+
+        if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
+        {
+            // 'degenerate' triangle check
+            valid[facei] = false;
+            hasInvalid = true;
+
+            if (verbose)
+            {
+                WarningInFunction
+                    << "triangle " << facei
+                    << " does not have three unique vertices:\n";
+                printTriangle(Warning, "    ", f, points());
+            }
+        }
+        else
+        {
+            // duplicate triangle check
+            const labelList& fEdges = faceEdges()[facei];
+
+            // Check if faceNeighbours use same points as this face.
+            // Note: discards normal information - sides of baffle are merged.
+
+            forAll(fEdges, fp)
+            {
+                const labelList& eFaces = edgeFaces()[fEdges[fp]];
+
+                forAll(eFaces, i)
+                {
+                    label neighbour = eFaces[i];
+
+                    if (neighbour > facei)
+                    {
+                        // lower numbered faces already checked
+                        const labelledTri& n = (*this)[neighbour];
+
+                        if
+                        (
+                            ((f[0] == n[0]) || (f[0] == n[1]) || (f[0] == n[2]))
+                         && ((f[1] == n[0]) || (f[1] == n[1]) || (f[1] == n[2]))
+                         && ((f[2] == n[0]) || (f[2] == n[1]) || (f[2] == n[2]))
+                        )
+                        {
+                            valid[facei] = false;
+                            hasInvalid = true;
+
+                            if (verbose)
+                            {
+                                WarningInFunction
+                                    << "triangles share the same vertices:\n"
+                                    << "    face 1 :" << facei << endl;
+                                printTriangle(Warning, "    ", f, points());
+
+                                Warning
+                                    << endl
+                                    << "    face 2 :"
+                                    << neighbour << endl;
+                                printTriangle(Warning, "    ", n, points());
+                            }
+
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    if (hasInvalid)
+    {
+        // Pack
+        label newFacei = 0;
+        forAll(*this, facei)
+        {
+            if (valid[facei])
+            {
+                const labelledTri& f = (*this)[facei];
+                (*this)[newFacei++] = f;
+            }
+        }
+
+        if (verbose)
+        {
+            WarningInFunction
+                << "Removing " << size() - newFacei
+                << " illegal faces." << endl;
+        }
+        (*this).setSize(newFacei);
+
+        // Topology can change because of renumbering
+        clearOut();
+    }
+}
+
+
+void Foam::triSurface::checkEdges(const bool verbose)
+{
+    const labelListList& eFaces = edgeFaces();
+
+    forAll(eFaces, edgeI)
+    {
+        const labelList& myFaces = eFaces[edgeI];
+
+        if (myFaces.empty())
+        {
+            FatalErrorInFunction
+                << "Edge " << edgeI << " with vertices " << edges()[edgeI]
+                << " has no edgeFaces"
+                << exit(FatalError);
+        }
+        else if (myFaces.size() > 2 && verbose)
+        {
+            WarningInFunction
+                << "Edge " << edgeI << " with vertices " << edges()[edgeI]
+                << " has more than 2 faces connected to it : " << myFaces
+                << endl;
+        }
+    }
+}
+
+
 void Foam::triSurface::cleanup(const bool verbose)
 {
     // Merge points (already done for STL, TRI)
@@ -827,8 +922,6 @@ void Foam::triSurface::cleanup(const bool verbose)
 }
 
 
-// Finds area, starting at facei, delimited by borderEdge. Marks all visited
-// faces (from face-edge-face walk) with currentZone.
 void Foam::triSurface::markZone
 (
     const boolList& borderEdge,
@@ -892,8 +985,6 @@ void Foam::triSurface::markZone
 }
 
 
-// Finds areas delimited by borderEdge (or 'real' edges).
-// Fills faceZone accordingly
 Foam::label Foam::triSurface::markZones
 (
     const boolList& borderEdge,
@@ -1040,64 +1131,177 @@ Foam::faceList Foam::triSurface::faces() const
 }
 
 
-Foam::scalar Foam::triSurface::pointNormalWeight
-(
-    const triFace& f,
-    const label pi,
-    const vector& fa,
-    const pointField& points
-) const
+Foam::tmp<Foam::scalarField> Foam::triSurface::curvature() const
 {
-    const 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(fa)/(magSqr(e1)*magSqr(e2) + vSmall);
-}
-
-
-Foam::tmp<Foam::vectorField> Foam::triSurface::pointNormals2() const
-{
-    // Weighted average of normals of faces attached to the vertex
-    // Weight = fA / (mag(e0)^2 * mag(e1)^2);
-
-    tmp<vectorField> tpointNormals
-    (
-        new vectorField(nPoints(), Zero)
-    );
-    vectorField& pointNormals = tpointNormals.ref();
+    // Curvature calculation is an implementation of the algorithm from:
+    // "Estimating Curvatures and their Derivatives on Triangle Meshes"
+    // by S. Rusinkiewicz
 
     const pointField& points = this->points();
-    const labelListList& pointFaces = this->pointFaces();
     const labelList& meshPoints = this->meshPoints();
+    const Map<label>& meshPointMap = this->meshPointMap();
 
-    forAll(pointFaces, pi)
+    const vectorField pointNormals
+    (
+        this->weightedPointNormals()
+    );
+
+    const triadField pointCoordSys
+    (
+        this->pointCoordSys(pointNormals)
+    );
+
+    tmp<scalarField> tcurvaturePointField(new scalarField(nPoints(), 0.0));
+    scalarField& curvaturePointField = tcurvaturePointField.ref();
+
+    List<symmTensor2D> pointFundamentalTensors(nPoints(), Zero);
+    scalarList accumulatedWeights(nPoints(), 0.0);
+
+    forAll(*this, fi)
     {
-        const labelList& pFaces = pointFaces[pi];
+        const triFace& f = operator[](fi);
+        const edgeList fEdges = f.edges();
 
-        forAll(pFaces, fi)
+        // Calculate the edge vectors and the normal differences
+        vectorField edgeVectors(f.size(), Zero);
+        vectorField normalDifferences(f.size(), Zero);
+
+        forAll(fEdges, feI)
         {
-            const label facei = pFaces[fi];
-            const triFace& f = operator[](facei);
+            const edge& e = fEdges[feI];
 
-            const vector fa = f.area(points);
-            const scalar weight =
-                pointNormalWeight(f, meshPoints[pi], fa, points);
-
-            pointNormals[pi] += weight*fa;
+            edgeVectors[feI] = e.vec(points);
+            normalDifferences[feI] =
+               pointNormals[meshPointMap[e[0]]]
+             - pointNormals[meshPointMap[e[1]]];
         }
 
-        pointNormals[pi] /= mag(pointNormals[pi]) + vSmall;
+        // Set up a local coordinate system for the face
+        const vector& e0 = edgeVectors[0];
+        const vector eN = f.area(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)
+        {
+            const scalar x = edgeVectors[i] & faceCoordSys[0];
+            const 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);
+
+            const scalar dndx = normalDifferences[i] & faceCoordSys[0];
+            const 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);
+        const 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
+            const tensor rotTensor = rotationTensor
+            (
+                ptCoordSys[2],
+                faceCoordSys[2]
+            );
+            const triad rotatedFaceCoordSys = rotTensor & tensor(faceCoordSys);
+
+            // Project the face curvature onto the point plane
+
+            const vector2D cmp1
+            (
+                ptCoordSys[0] & rotatedFaceCoordSys[0],
+                ptCoordSys[0] & rotatedFaceCoordSys[1]
+            );
+            const vector2D cmp2
+            (
+                ptCoordSys[1] & rotatedFaceCoordSys[0],
+                ptCoordSys[1] & rotatedFaceCoordSys[1]
+            );
+
+            const tensor2D projTensor
+            (
+                cmp1,
+                cmp2
+            );
+
+            const symmTensor2D projectedFundamentalTensor
+            (
+                projTensor.x() & (secondFundamentalTensor & projTensor.x()),
+                projTensor.x() & (secondFundamentalTensor & projTensor.y()),
+                projTensor.y() & (secondFundamentalTensor & projTensor.y())
+            );
+
+            // Calculate weight
+            const scalar weight = pointNormalWeight
+            (
+                f,
+                meshPoints[patchPointIndex],
+                f.area(points),
+                points
+            );
+
+            // Sum contribution of face to this point
+            pointFundamentalTensors[patchPointIndex] +=
+                weight*projectedFundamentalTensor;
+
+            accumulatedWeights[patchPointIndex] += weight;
+        }
     }
 
-    return tpointNormals;
+    forAll(curvaturePointField, pi)
+    {
+        pointFundamentalTensors[pi] /= (accumulatedWeights[pi] + small);
+
+        const vector2D principalCurvatures =
+            eigenValues(pointFundamentalTensors[pi]);
+
+        //scalar curvature =
+        //    (principalCurvatures[0] + principalCurvatures[1])/2;
+        const scalar curvature = max
+        (
+            mag(principalCurvatures[0]),
+            mag(principalCurvatures[1])
+        );
+        //scalar curvature = principalCurvatures[0]*principalCurvatures[1];
+
+        curvaturePointField[meshPoints[pi]] = curvature;
+    }
+
+    return tcurvaturePointField;
 }
 
 
diff --git a/src/triSurface/triSurface/triSurface.H b/src/triSurface/triSurface/triSurface.H
index 6699e84b9..d33ef6a7f 100644
--- a/src/triSurface/triSurface/triSurface.H
+++ b/src/triSurface/triSurface/triSurface.H
@@ -42,7 +42,7 @@ SourceFiles
 #include "geometricSurfacePatchList.H"
 #include "surfacePatchList.H"
 #include "triFaceList.H"
-#include "triSurfaceFieldsFwd.H"
+#include "triadField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -125,6 +125,21 @@ class triSurface
             const bool verbose = false
         );
 
+        scalar pointNormalWeight
+        (
+            const triFace& f,
+            const label pi,
+            const vector& fa,
+            const pointField& points
+        ) const;
+
+        //- Return the surface point normals
+        tmp<vectorField> weightedPointNormals() const;
+
+        //- Return the curvature of surface at the points
+        tmp<triadField> pointCoordSys(const vectorField& pointNormals) const;
+
+
         //- Read in Foam format
         bool read(Istream&);
 
@@ -397,22 +412,8 @@ public:
 
         // Analysis
 
-            scalar pointNormalWeight
-            (
-                const triFace& f,
-                const label pi,
-                const vector& fa,
-                const pointField& points
-            ) const;
-
-            //- Return the surface point normals
-            tmp<vectorField> pointNormals2() const;
-
             //- Return the curvature of surface at the points
-            tmp<triSurfacePointScalarField> pointCoordSys() const;
-
-            //- Return the curvature of surface at the points
-            tmp<triSurfacePointScalarField> curvature() const;
+            tmp<scalarField> curvature() const;
 
 
         // Write