foamyHexMesh: Added cell sizing based on local surface closeness

First run the surfaceFeatureExtract with the "closeness" option enabled in the
surfaceFeatureExtractDict to extract the surface closeness point field

    // Out put the closeness of surface elements to other surface elements.
    closeness               yes;

Then enable cell sizing based on local surface closeness by specifying the
"internalCloseness" options in the foamyHexMeshDict e.g.

motionControl
{
    defaultCellSize             4;

    minimumCellSizeCoeff        0.1;
    maxSmoothingIterations      100;
    maxRefinementIterations     2;

    shapeControlFunctions
    {
        geometry
        {
            type                        searchableSurfaceControl;
            priority                    1;
            mode                        inside;

            surfaceCellSizeFunction     nonUniformField;

            cellSizeCalculationType     automatic;

            curvature                   false;
            curvatureFile               dummy;
            featureProximity            false;
            featureProximityFile        dummy;
            internalCloseness           true;
            internalClosenessFile       geometry.internalPointCloseness;
            internalClosenessCellSizeCoeff 25;
            curvatureCellSizeCoeff      0;
            maximumCellSizeCoeff        1;
            cellSizeFunction            uniform;
        }
    }
}

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,6 +115,10 @@ Foam::automatic::automatic
     featureProximityFile_(coeffsDict_.lookup("featureProximityFile")),
     readInternalCloseness_(Switch(coeffsDict_.lookup("internalCloseness"))),
     internalClosenessFile_(coeffsDict_.lookup("internalClosenessFile")),
+    internalClosenessCellSizeCoeff_
+    (
+        readScalar(coeffsDict_.lookup("internalClosenessCellSizeCoeff"))
+    ),
     curvatureCellSizeCoeff_
     (
         readScalar(coeffsDict_.lookup("curvatureCellSizeCoeff"))
@@ -203,7 +207,7 @@ Foam::tmp<Foam::triSurfacePointScalarField> Foam::automatic::load()
         Info<< indent
             << "Reading internal closeness: " << internalClosenessFile_ << endl;
 
-        triSurfaceScalarField internalCloseness
+        triSurfacePointScalarField internalClosenessPointField
         (
             IOobject
             (
@@ -219,17 +223,13 @@ Foam::tmp<Foam::triSurfacePointScalarField> Foam::automatic::load()
             true
         );
 
-        scalarField internalClosenessPointField
-        (
-            patchInterpolate.faceToPointInterpolate(internalCloseness)
-        );
-
         forAll(pointCellSize, pI)
         {
             pointCellSize[pI] =
                 min
                 (
-                    internalClosenessPointField[meshPointMap[pI]],
+                    (1.0/internalClosenessCellSizeCoeff_)
+                   *internalClosenessPointField[meshPointMap[pI]],
                     pointCellSize[pI]
                 );
         }

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellSizeControlSurfaces/surfaceCellSizeFunction/cellSizeCalculationType/automatic/automatic.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -74,6 +74,10 @@ private:
         const Switch readInternalCloseness_;
         const word internalClosenessFile_;
 
+        //- The internalCloseness values are multiplied by the inverse
+        //  of this value to get the cell size
+        const scalar internalClosenessCellSizeCoeff_;
+
         //- The curvature values are multiplied by the inverse of this value to
         //  get the cell size
         const scalar curvatureCellSizeCoeff_;

applications/utilities/surface/surfaceFeatureExtract/Make/files

@@ -1,3 +1,6 @@
+surfaceExtractCloseness.C
+surfaceExtractPointCloseness.C
+surfaceFeatureExtractUtilities.C
 surfaceFeatureExtract.C
 
 EXE = $(FOAM_APPBIN)/surfaceFeatureExtract

applications/utilities/surface/surfaceFeatureExtract/surfaceExtractCloseness.C

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeatureExtract.H"
+#include "Time.H"
+#include "triSurfaceMesh.H"
+#include "vtkSurfaceWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::extractCloseness
+(
+    const fileName &sFeatFileName,
+    const Time& runTime,
+    const triSurface &surf,
+    const bool writeVTK
+)
+{
+    // Searchable triSurface
+    const triSurfaceMesh searchSurf
+    (
+        IOobject
+        (
+            sFeatFileName + ".closeness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf
+    );
+
+
+    // Prepare start and end points for intersection tests
+
+    const vectorField& normals = searchSurf.faceNormals();
+
+    const scalar span = searchSurf.bounds().mag();
+
+    const pointField start(searchSurf.faceCentres() - span*normals);
+    const pointField end(searchSurf.faceCentres() + span*normals);
+    const pointField& faceCentres = searchSurf.faceCentres();
+
+    List<List<pointIndexHit>> allHitinfo;
+
+    // Find all intersections (in order)
+    searchSurf.findLineAll(start, end, allHitinfo);
+
+    scalarField internalCloseness(start.size(), GREAT);
+    scalarField externalCloseness(start.size(), GREAT);
+
+    forAll(allHitinfo, fi)
+    {
+        const List<pointIndexHit>& hitInfo = allHitinfo[fi];
+
+        processHit
+        (
+            internalCloseness[fi],
+            externalCloseness[fi],
+            fi,
+            surf,
+            start[fi],
+            faceCentres[fi],
+            end[fi],
+            normals[fi],
+            normals,
+            hitInfo
+        );
+    }
+
+    triSurfaceScalarField internalClosenessField
+    (
+        IOobject
+        (
+            sFeatFileName + ".internalCloseness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf,
+        dimLength,
+        internalCloseness
+    );
+
+    internalClosenessField.write();
+
+    triSurfaceScalarField externalClosenessField
+    (
+        IOobject
+        (
+            sFeatFileName + ".externalCloseness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf,
+        dimLength,
+        externalCloseness
+    );
+
+    externalClosenessField.write();
+
+    if (writeVTK)
+    {
+        const faceList faces(surf.faces());
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constantPath()/"triSurface",// outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "internalCloseness",                // fieldName
+            internalCloseness,
+            false,                              // isNodeValues
+            true                                // verbose
+        );
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constantPath()/"triSurface",// outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "externalCloseness",                // fieldName
+            externalCloseness,
+            false,                              // isNodeValues
+            true                                // verbose
+        );
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/surfaceExtractPointCloseness.C

@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeatureExtract.H"
+#include "Time.H"
+#include "triSurfaceMesh.H"
+#include "vtkSurfaceWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void Foam::processHit
+(
+    scalar& internalCloseness,
+    scalar& externalCloseness,
+    const label fi,
+    const triSurface& surf,
+    const point& start,
+    const point& p,
+    const point& end,
+    const vector& normal,
+    const vectorField& normals,
+    const List<pointIndexHit>& hitInfo
+)
+{
+    if (hitInfo.size() < 1)
+    {
+        drawHitProblem(fi, surf, start, p, end, hitInfo);
+    }
+    else if (hitInfo.size() == 1)
+    {
+        if (!hitInfo[0].hit())
+        {
+        }
+        else if (hitInfo[0].index() != fi)
+        {
+            drawHitProblem(fi, surf, start, p, end, hitInfo);
+        }
+    }
+    else
+    {
+        label ownHiti = -1;
+
+        forAll(hitInfo, hI)
+        {
+            // Find the hit on the triangle that launched the ray
+
+            if (hitInfo[hI].index() == fi)
+            {
+                ownHiti = hI;
+                break;
+            }
+        }
+
+        if (ownHiti < 0)
+        {
+            drawHitProblem(fi, surf, start, p, end, hitInfo);
+        }
+        else if (ownHiti == 0)
+        {
+            // There are no internal hits, the first hit is the
+            // closest external hit
+
+            if
+            (
+                (normal & normals[hitInfo[ownHiti + 1].index()])
+              < externalToleranceCosAngle
+            )
+            {
+                externalCloseness = min
+                (
+                    externalCloseness,
+                    mag(p - hitInfo[ownHiti + 1].hitPoint())
+                );
+            }
+        }
+        else if (ownHiti == hitInfo.size() - 1)
+        {
+            // There are no external hits, the last but one hit is
+            // the closest internal hit
+
+            if
+            (
+                (normal & normals[hitInfo[ownHiti - 1].index()])
+              < internalToleranceCosAngle
+            )
+            {
+                internalCloseness = min
+                (
+                    internalCloseness,
+                    mag(p - hitInfo[ownHiti - 1].hitPoint())
+                );
+            }
+        }
+        else
+        {
+            if
+            (
+                (normal & normals[hitInfo[ownHiti + 1].index()])
+              < externalToleranceCosAngle
+            )
+            {
+                externalCloseness = min
+                (
+                    externalCloseness,
+                    mag(p - hitInfo[ownHiti + 1].hitPoint())
+                );
+            }
+
+            if
+            (
+                (normal & normals[hitInfo[ownHiti - 1].index()])
+              < internalToleranceCosAngle
+            )
+            {
+                internalCloseness = min
+                (
+                    internalCloseness,
+                    mag(p - hitInfo[ownHiti - 1].hitPoint())
+                );
+            }
+        }
+    }
+}
+
+
+void Foam::extractPointCloseness
+(
+    const fileName &sFeatFileName,
+    const Time& runTime,
+    const triSurface &surf,
+    const bool writeVTK
+)
+{
+    // Searchable triSurface
+    const triSurfaceMesh searchSurf
+    (
+        IOobject
+        (
+            sFeatFileName + ".closeness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf
+    );
+
+
+    // Prepare start and end points for intersection tests
+
+    const pointField& points = searchSurf.points();
+    const labelList& meshPoints = searchSurf.meshPoints();
+    const pointField& faceCentres = searchSurf.faceCentres();
+    const vectorField& normals = searchSurf.faceNormals();
+    const labelListList& pointFaces = searchSurf.pointFaces();
+
+    const scalar span = searchSurf.bounds().mag();
+
+    label nPointFaces = 0;
+    forAll(pointFaces, pfi)
+    {
+        nPointFaces += pointFaces[pfi].size();
+    }
+
+    pointField facePoints(nPointFaces);
+    pointField start(nPointFaces);
+    pointField end(nPointFaces);
+
+    label i = 0;
+    forAll(points, pi)
+    {
+        forAll(pointFaces[pi], pfi)
+        {
+            const label fi = pointFaces[pi][pfi];
+
+            facePoints[i] = (0.9*points[meshPoints[pi]] + 0.1*faceCentres[fi]);
+            const vector& n = normals[fi];
+
+            start[i] = facePoints[i] - span*n;
+            end[i] = facePoints[i] + span*n;
+
+            i++;
+        }
+    }
+
+    List<List<pointIndexHit>> allHitinfo;
+
+    // Find all intersections (in order)
+    searchSurf.findLineAll(start, end, allHitinfo);
+
+    scalarField internalCloseness(points.size(), GREAT);
+    scalarField externalCloseness(points.size(), GREAT);
+
+    i = 0;
+    forAll(points, pi)
+    {
+        forAll(pointFaces[pi], pfi)
+        {
+            const label fi = pointFaces[pi][pfi];
+            const List<pointIndexHit>& hitInfo = allHitinfo[i];
+
+            processHit
+            (
+                internalCloseness[pi],
+                externalCloseness[pi],
+                fi,
+                surf,
+                start[i],
+                facePoints[i],
+                end[i],
+                normals[fi],
+                normals,
+                hitInfo
+            );
+
+            i++;
+        }
+    }
+
+    triSurfacePointScalarField internalClosenessPointField
+    (
+        IOobject
+        (
+            sFeatFileName + ".internalPointCloseness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf,
+        dimLength,
+        internalCloseness
+    );
+
+    internalClosenessPointField.write();
+
+    triSurfacePointScalarField externalClosenessPointField
+    (
+        IOobject
+        (
+            sFeatFileName + ".externalPointCloseness",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        surf,
+        dimLength,
+        externalCloseness
+    );
+
+    externalClosenessPointField.write();
+
+    if (writeVTK)
+    {
+        const faceList faces(surf.faces());
+        const Map<label>& meshPointMap = surf.meshPointMap();
+
+        forAll(meshPointMap, pi)
+        {
+            internalCloseness[pi] =
+                internalClosenessPointField[meshPointMap[pi]];
+
+            externalCloseness[pi] =
+                externalClosenessPointField[meshPointMap[pi]];
+        }
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constantPath()/"triSurface",// outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "internalPointCloseness",           // fieldName
+            internalCloseness,
+            true,                               // isNodeValues
+            true                                // verbose
+        );
+
+        vtkSurfaceWriter().write
+        (
+            runTime.constantPath()/"triSurface",// outputDir
+            sFeatFileName,                      // surfaceName
+            surf.points(),
+            faces,
+            "externalPointCloseness",           // fieldName
+            externalCloseness,
+            true,                               // isNodeValues
+            true                                // verbose
+        );
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.H

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    surfaceFeatureExtract
+
+Description
+    Utility functions for surfaceFeatureExtract
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeatures.H"
+#include "extendedFeatureEdgeMesh.H"
+#include "triSurfaceFields.H"
+#include "plane.H"
+#include "triadField.H"
+
+namespace Foam
+{
+    extern const scalar internalAngleTolerance;
+    extern const scalar internalToleranceCosAngle;
+
+    extern const scalar externalAngleTolerance;
+    extern const scalar externalToleranceCosAngle;
+
+    scalar calcVertexNormalWeight
+    (
+        const triFace& f,
+        const label pI,
+        const vector& fN,
+        const pointField& points
+    );
+
+    point randomPointInPlane(const plane& p);
+
+    triadField calcVertexCoordSys
+    (
+        const triSurface& surf,
+        const vectorField& pointNormals
+    );
+
+    vectorField calcVertexNormals(const triSurface& surf);
+
+    triSurfacePointScalarField calcCurvature
+    (
+        const word& name,
+        const Time& runTime,
+        const triSurface& surf,
+        const vectorField& pointNormals,
+        const triadField& pointCoordSys
+    );
+
+    bool edgesConnected(const edge& e1, const edge& e2);
+
+    scalar calcProximityOfFeaturePoints
+    (
+        const List<pointIndexHit>& hitList,
+        const scalar defaultCellSize
+    );
+
+    scalar calcProximityOfFeatureEdges
+    (
+        const extendedFeatureEdgeMesh& efem,
+        const List<pointIndexHit>& hitList,
+        const scalar defaultCellSize
+    );
+
+    void dumpBox(const treeBoundBox& bb, const fileName& fName);
+
+    //- Deletes all edges inside/outside bounding box from set.
+    void deleteBox
+    (
+        const triSurface& surf,
+        const treeBoundBox& bb,
+        const bool removeInside,
+        List<surfaceFeatures::edgeStatus>& edgeStat
+    );
+
+    bool onLine(const point& p, const linePointRef& line);
+
+    //- Deletes all edges inside/outside bounding box from set.
+    void deleteEdges
+    (
+        const triSurface& surf,
+        const plane& cutPlane,
+        List<surfaceFeatures::edgeStatus>& edgeStat
+    );
+
+    void processHit
+    (
+        scalar& internalCloseness,
+        scalar& externalCloseness,
+        const label fi,
+        const triSurface& surf,
+        const point& start,
+        const point& p,
+        const point& end,
+        const vector& normal,
+        const vectorField& normals,
+        const List<pointIndexHit>& hitInfo
+    );
+
+    void drawHitProblem
+    (
+        const label fi,
+        const triSurface& surf,
+        const point& start,
+        const point& p,
+        const point& end,
+        const List<pointIndexHit>& hitInfo
+    );
+
+    //- Unmark non-manifold edges if individual triangles are not features
+    void unmarkBaffles
+    (
+        const triSurface& surf,
+        const scalar includedAngle,
+        List<surfaceFeatures::edgeStatus>& edgeStat
+    );
+
+    //- Divide into multiple normal bins
+    //  - return REGION if != 2 normals
+    //  - return REGION if 2 normals that make feature angle
+    //  - otherwise return NONE and set normals,bins
+    surfaceFeatures::edgeStatus checkFlatRegionEdge
+    (
+        const triSurface& surf,
+        const scalar tol,
+        const scalar includedAngle,
+        const label edgeI
+    );
+
+    void extractCloseness
+    (
+        const fileName &sFeatFileName,
+        const Time& runTime,
+        const triSurface &surf,
+        const bool writeVTK
+    );
+
+    void extractPointCloseness
+    (
+        const fileName &sFeatFileName,
+        const Time& runTime,
+        const triSurface &surf,
+        const bool writeVTK
+    );
+
+    void writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os);
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtractUtilities.C

@@ -0,0 +1,938 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    surfaceFeatureExtract
+
+Description
+    Extracts and writes surface features to file. All but the basic feature
+    extraction is WIP.
+
+    Curvature calculation is an implementation of the algorithm from:
+
+      "Estimating Curvatures and their Derivatives on Triangle Meshes"
+      by S. Rusinkiewicz
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceFeatureExtract.H"
+#include "Time.H"
+#include "meshTools.H"
+#include "tensor2D.H"
+#include "symmTensor2D.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const Foam::scalar Foam::internalAngleTolerance(45);
+const Foam::scalar Foam::internalToleranceCosAngle
+(
+    cos(degToRad(180 - internalAngleTolerance))
+);
+
+const Foam::scalar Foam::externalAngleTolerance(10);
+const Foam::scalar Foam::externalToleranceCosAngle
+(
+    cos(degToRad(180 - externalAngleTolerance))
+);
+
+
+Foam::scalar Foam::calcVertexNormalWeight
+(
+    const triFace& f,
+    const label pI,
+    const vector& fN,
+    const pointField& points
+)
+{
+    label index = findIndex(f, pI);
+
+    if (index == -1)
+    {
+        FatalErrorInFunction
+            << "Point not in face" << abort(FatalError);
+    }
+
+    const vector e1 = points[f[index]] - points[f[f.fcIndex(index)]];
+    const vector e2 = points[f[index]] - points[f[f.rcIndex(index)]];
+
+    return mag(fN)/(magSqr(e1)*magSqr(e2) + VSMALL);
+}
+
+
+Foam::point Foam::randomPointInPlane(const plane& p)
+{
+    // Perturb base point
+    const point& refPt = p.refPoint();
+
+    // ax + by + cz + d = 0
+    const FixedList<scalar, 4>& planeCoeffs = p.planeCoeffs();
+
+    const scalar perturbX = refPt.x() + 1e-3;
+    const scalar perturbY = refPt.y() + 1e-3;
+    const scalar perturbZ = refPt.z() + 1e-3;
+
+    if (mag(planeCoeffs[2]) < SMALL)
+    {
+        if (mag(planeCoeffs[1]) < SMALL)
+        {
+            const scalar x =
+                -1.0
+                *(
+                     planeCoeffs[3]
+                   + planeCoeffs[1]*perturbY
+                   + planeCoeffs[2]*perturbZ
+                 )/planeCoeffs[0];
+
+            return point
+            (
+                x,
+                perturbY,
+                perturbZ
+            );
+        }
+
+        const scalar y =
+            -1.0
+            *(
+                 planeCoeffs[3]
+               + planeCoeffs[0]*perturbX
+               + planeCoeffs[2]*perturbZ
+             )/planeCoeffs[1];
+
+        return point
+        (
+            perturbX,
+            y,
+            perturbZ
+        );
+    }
+    else
+    {
+        const scalar z =
+            -1.0
+            *(
+                 planeCoeffs[3]
+               + planeCoeffs[0]*perturbX
+               + planeCoeffs[1]*perturbY
+             )/planeCoeffs[2];
+
+        return point
+        (
+            perturbX,
+            perturbY,
+            z
+        );
+    }
+}
+
+
+Foam::triadField Foam::calcVertexCoordSys
+(
+    const triSurface& surf,
+    const vectorField& pointNormals
+)
+{
+    const pointField& points = surf.points();
+    const Map<label>& meshPointMap = surf.meshPointMap();
+
+    triadField pointCoordSys(points.size());
+
+    forAll(points, pI)
+    {
+        const point& pt = points[pI];
+        const vector& normal = pointNormals[meshPointMap[pI]];
+
+        if (mag(normal) < SMALL)
+        {
+            pointCoordSys[meshPointMap[pI]] = triad::unset;
+            continue;
+        }
+
+        plane p(pt, normal);
+
+        // Pick random point in plane
+        vector dir1 = pt - randomPointInPlane(p);
+        dir1 /= mag(dir1);
+
+        vector dir2 = dir1 ^ normal;
+        dir2 /= mag(dir2);
+
+        pointCoordSys[meshPointMap[pI]] = triad(dir1, dir2, normal);
+    }
+
+    return pointCoordSys;
+}
+
+
+Foam::vectorField Foam::calcVertexNormals(const triSurface& surf)
+{
+    // Weighted average of normals of faces attached to the vertex
+    // Weight = fA / (mag(e0)^2 * mag(e1)^2);
+
+    Info<< "Calculating vertex normals" << endl;
+
+    vectorField pointNormals(surf.nPoints(), Zero);
+
+    const pointField& points = surf.points();
+    const labelListList& pointFaces = surf.pointFaces();
+    const labelList& meshPoints = surf.meshPoints();
+
+    forAll(pointFaces, pI)
+    {
+        const labelList& pFaces = pointFaces[pI];
+
+        forAll(pFaces, fi)
+        {
+            const label facei = pFaces[fi];
+            const triFace& f = surf[facei];
+
+            vector fN = f.normal(points);
+
+            scalar weight = calcVertexNormalWeight
+            (
+                f,
+                meshPoints[pI],
+                fN,
+                points
+            );
+
+            pointNormals[pI] += weight*fN;
+        }
+
+        pointNormals[pI] /= mag(pointNormals[pI]) + VSMALL;
+    }
+
+    return pointNormals;
+}
+
+
+Foam::triSurfacePointScalarField Foam::calcCurvature
+(
+    const word& name,
+    const Time& runTime,
+    const triSurface& surf,
+    const vectorField& pointNormals,
+    const triadField& pointCoordSys
+)
+{
+    Info<< "Calculating face curvature" << endl;
+
+    const pointField& points = surf.points();
+    const labelList& meshPoints = surf.meshPoints();
+    const Map<label>& meshPointMap = surf.meshPointMap();
+
+    triSurfacePointScalarField curvaturePointField
+    (
+        IOobject
+        (
+            name + ".curvature",
+            runTime.constant(),
+            "triSurface",
+            runTime,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        surf,
+        dimLength,
+        scalarField(points.size(), 0.0)
+    );
+
+    List<symmTensor2D> pointFundamentalTensors
+    (
+        points.size(),
+        symmTensor2D::zero
+    );
+
+    scalarList accumulatedWeights(points.size(), 0.0);
+
+    forAll(surf, fi)
+    {
+        const triFace& f = surf[fi];
+        const edgeList fEdges = f.edges();
+
+        // Calculate the edge vectors and the normal differences
+        vectorField edgeVectors(f.size(), Zero);
+        vectorField normalDifferences(f.size(), Zero);
+
+        forAll(fEdges, feI)
+        {
+            const edge& e = fEdges[feI];
+
+            edgeVectors[feI] = e.vec(points);
+            normalDifferences[feI] =
+               pointNormals[meshPointMap[e[0]]]
+             - pointNormals[meshPointMap[e[1]]];
+        }
+
+        // Set up a local coordinate system for the face
+        const vector& e0 = edgeVectors[0];
+        const vector eN = f.normal(points);
+        const vector e1 = (e0 ^ eN);
+
+        if (magSqr(eN) < ROOTVSMALL)
+        {
+            continue;
+        }
+
+        triad faceCoordSys(e0, e1, eN);
+        faceCoordSys.normalize();
+
+        // Construct the matrix to solve
+        scalarSymmetricSquareMatrix T(3, 0);
+        scalarDiagonalMatrix Z(3, 0);
+
+        // Least Squares
+        for (label i = 0; i < 3; ++i)
+        {
+            scalar x = edgeVectors[i] & faceCoordSys[0];
+            scalar y = edgeVectors[i] & faceCoordSys[1];
+
+            T(0, 0) += sqr(x);
+            T(1, 0) += x*y;
+            T(1, 1) += sqr(x) + sqr(y);
+            T(2, 1) += x*y;
+            T(2, 2) += sqr(y);
+
+            scalar dndx = normalDifferences[i] & faceCoordSys[0];
+            scalar dndy = normalDifferences[i] & faceCoordSys[1];
+
+            Z[0] += dndx*x;
+            Z[1] += dndx*y + dndy*x;
+            Z[2] += dndy*y;
+        }
+
+        // Perform Cholesky decomposition and back substitution.
+        // Decomposed matrix is in T and solution is in Z.
+        LUsolve(T, Z);
+        symmTensor2D secondFundamentalTensor(Z[0], Z[1], Z[2]);
+
+        // Loop over the face points adding the contribution of the face
+        // curvature to the points.
+        forAll(f, fpI)
+        {
+            const label patchPointIndex = meshPointMap[f[fpI]];
+
+            const triad& ptCoordSys = pointCoordSys[patchPointIndex];
+
+            if (!ptCoordSys.set())
+            {
+                continue;
+            }
+
+            // Rotate faceCoordSys to ptCoordSys
+            tensor rotTensor = rotationTensor(ptCoordSys[2], faceCoordSys[2]);
+            triad rotatedFaceCoordSys = rotTensor & tensor(faceCoordSys);
+
+            // Project the face curvature onto the point plane
+
+            vector2D cmp1
+            (
+                ptCoordSys[0] & rotatedFaceCoordSys[0],
+                ptCoordSys[0] & rotatedFaceCoordSys[1]
+            );
+            vector2D cmp2
+            (
+                ptCoordSys[1] & rotatedFaceCoordSys[0],
+                ptCoordSys[1] & rotatedFaceCoordSys[1]
+            );
+
+            tensor2D projTensor
+            (
+                cmp1,
+                cmp2
+            );
+
+            symmTensor2D projectedFundamentalTensor
+            (
+                projTensor.x() & (secondFundamentalTensor & projTensor.x()),
+                projTensor.x() & (secondFundamentalTensor & projTensor.y()),
+                projTensor.y() & (secondFundamentalTensor & projTensor.y())
+            );
+
+            // Calculate weight
+            // TODO: Voronoi area weighting
+            scalar weight = calcVertexNormalWeight
+            (
+                f,
+                meshPoints[patchPointIndex],
+                f.normal(points),
+                points
+            );
+
+            // Sum contribution of face to this point
+            pointFundamentalTensors[patchPointIndex] +=
+                weight*projectedFundamentalTensor;
+
+            accumulatedWeights[patchPointIndex] += weight;
+        }
+
+        if (false)
+        {
+            Info<< "Points = " << points[f[0]] << " "
+                << points[f[1]] << " "
+                << points[f[2]] << endl;
+            Info<< "edgeVecs = " << edgeVectors[0] << " "
+                << edgeVectors[1] << " "
+                << edgeVectors[2] << endl;
+            Info<< "normDiff = " << normalDifferences[0] << " "
+                << normalDifferences[1] << " "
+                << normalDifferences[2] << endl;
+            Info<< "faceCoordSys = " << faceCoordSys << endl;
+            Info<< "T = " << T << endl;
+            Info<< "Z = " << Z << endl;
+        }
+    }
+
+    forAll(curvaturePointField, pI)
+    {
+        pointFundamentalTensors[pI] /= (accumulatedWeights[pI] + SMALL);
+
+        vector2D principalCurvatures = eigenValues(pointFundamentalTensors[pI]);
+
+        //scalar curvature =
+        //    (principalCurvatures[0] + principalCurvatures[1])/2;
+        scalar curvature = max
+        (
+            mag(principalCurvatures[0]),
+            mag(principalCurvatures[1])
+        );
+        //scalar curvature = principalCurvatures[0]*principalCurvatures[1];
+
+        curvaturePointField[meshPoints[pI]] = curvature;
+    }
+
+    return curvaturePointField;
+}
+
+
+bool Foam::edgesConnected(const edge& e1, const edge& e2)
+{
+    if
+    (
+        e1.start() == e2.start()
+     || e1.start() == e2.end()
+     || e1.end() == e2.start()
+     || e1.end() == e2.end()
+    )
+    {
+        return true;
+    }
+
+    return false;
+}
+
+
+Foam::scalar Foam::calcProximityOfFeaturePoints
+(
+    const List<pointIndexHit>& hitList,
+    const scalar defaultCellSize
+)
+{
+    scalar minDist = defaultCellSize;
+
+    for
+    (
+        label hI1 = 0;
+        hI1 < hitList.size() - 1;
+        ++hI1
+    )
+    {
+        const pointIndexHit& pHit1 = hitList[hI1];
+
+        if (pHit1.hit())
+        {
+            for
+            (
+                label hI2 = hI1 + 1;
+                hI2 < hitList.size();
+                ++hI2
+            )
+            {
+                const pointIndexHit& pHit2 = hitList[hI2];
+
+                if (pHit2.hit())
+                {
+                    scalar curDist = mag(pHit1.hitPoint() - pHit2.hitPoint());
+
+                    minDist = min(curDist, minDist);
+                }
+            }
+        }
+    }
+
+    return minDist;
+}
+
+
+Foam::scalar Foam::calcProximityOfFeatureEdges
+(
+    const extendedFeatureEdgeMesh& efem,
+    const List<pointIndexHit>& hitList,
+    const scalar defaultCellSize
+)
+{
+    scalar minDist = defaultCellSize;
+
+    for
+    (
+        label hI1 = 0;
+        hI1 < hitList.size() - 1;
+        ++hI1
+    )
+    {
+        const pointIndexHit& pHit1 = hitList[hI1];
+
+        if (pHit1.hit())
+        {
+            const edge& e1 = efem.edges()[pHit1.index()];
+
+            for
+            (
+                label hI2 = hI1 + 1;
+                hI2 < hitList.size();
+                ++hI2
+            )
+            {
+                const pointIndexHit& pHit2 = hitList[hI2];
+
+                if (pHit2.hit())
+                {
+                    const edge& e2 = efem.edges()[pHit2.index()];
+
+                    // Don't refine if the edges are connected to each other
+                    if (!edgesConnected(e1, e2))
+                    {
+                        scalar curDist =
+                            mag(pHit1.hitPoint() - pHit2.hitPoint());
+
+                        minDist = min(curDist, minDist);
+                    }
+                }
+            }
+        }
+    }
+
+    return minDist;
+}
+
+
+void Foam::dumpBox(const treeBoundBox& bb, const fileName& fName)
+{
+    OFstream str(fName);
+
+    Info<< "Dumping bounding box " << bb << " as lines to obj file "
+        << str.name() << endl;
+
+
+    pointField boxPoints(bb.points());
+
+    forAll(boxPoints, i)
+    {
+        meshTools::writeOBJ(str, boxPoints[i]);
+    }
+
+    forAll(treeBoundBox::edges, i)
+    {
+        const edge& e = treeBoundBox::edges[i];
+
+        str<< "l " << e[0]+1 <<  ' ' << e[1]+1 << nl;
+    }
+}
+
+
+void Foam::deleteBox
+(
+    const triSurface& surf,
+    const treeBoundBox& bb,
+    const bool removeInside,
+    List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+    forAll(edgeStat, edgeI)
+    {
+        const point eMid = surf.edges()[edgeI].centre(surf.localPoints());
+
+        if (removeInside ? bb.contains(eMid) : !bb.contains(eMid))
+        {
+            edgeStat[edgeI] = surfaceFeatures::NONE;
+        }
+    }
+}
+
+
+bool Foam::onLine(const point& p, const linePointRef& line)
+{
+    const point& a = line.start();
+    const point& b = line.end();
+
+    if
+    (
+        ( p.x() < min(a.x(), b.x()) || p.x() > max(a.x(), b.x()) )
+     || ( p.y() < min(a.y(), b.y()) || p.y() > max(a.y(), b.y()) )
+     || ( p.z() < min(a.z(), b.z()) || p.z() > max(a.z(), b.z()) )
+    )
+    {
+        return false;
+    }
+
+    return true;
+}
+
+
+void Foam::deleteEdges
+(
+    const triSurface& surf,
+    const plane& cutPlane,
+    List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+    const pointField& points = surf.points();
+    const labelList& meshPoints = surf.meshPoints();
+
+    forAll(edgeStat, edgeI)
+    {
+        const edge& e = surf.edges()[edgeI];
+        const point& p0 = points[meshPoints[e.start()]];
+        const point& p1 = points[meshPoints[e.end()]];
+        const linePointRef line(p0, p1);
+
+        // If edge does not intersect the plane, delete.
+        scalar intersect = cutPlane.lineIntersect(line);
+
+        point featPoint = intersect * (p1 - p0) + p0;
+
+        if (!onLine(featPoint, line))
+        {
+            edgeStat[edgeI] = surfaceFeatures::NONE;
+        }
+    }
+}
+
+
+void Foam::drawHitProblem
+(
+    const label fi,
+    const triSurface& surf,
+    const point& start,
+    const point& p,
+    const point& end,
+    const List<pointIndexHit>& hitInfo
+)
+{
+    Info<< nl << "# findLineAll did not hit its own face."
+        << nl << "# fi " << fi
+        << nl << "# start " << start
+        << nl << "# point " << p
+        << nl << "# end " << end
+        << nl << "# hitInfo " << hitInfo
+        << endl;
+
+    meshTools::writeOBJ(Info, start);
+    meshTools::writeOBJ(Info, p);
+    meshTools::writeOBJ(Info, end);
+
+    Info<< "l 1 2 3" << endl;
+
+    meshTools::writeOBJ(Info, surf.points()[surf[fi][0]]);
+    meshTools::writeOBJ(Info, surf.points()[surf[fi][1]]);
+    meshTools::writeOBJ(Info, surf.points()[surf[fi][2]]);
+
+    Info<< "f 4 5 6" << endl;
+
+    forAll(hitInfo, hI)
+    {
+        label hFI = hitInfo[hI].index();
+
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][0]]);
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][1]]);
+        meshTools::writeOBJ(Info, surf.points()[surf[hFI][2]]);
+
+        Info<< "f "
+            << 3*hI + 7 << " "
+            << 3*hI + 8 << " "
+            << 3*hI + 9
+            << endl;
+    }
+}
+
+
+void Foam::unmarkBaffles
+(
+    const triSurface& surf,
+    const scalar includedAngle,
+    List<surfaceFeatures::edgeStatus>& edgeStat
+)
+{
+    scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
+
+    const labelListList& edgeFaces = surf.edgeFaces();
+
+    forAll(edgeFaces, edgeI)
+    {
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        if (eFaces.size() > 2)
+        {
+            label i0 = eFaces[0];
+            //const labelledTri& f0 = surf[i0];
+            const Foam::vector& n0 = surf.faceNormals()[i0];
+
+            //Pout<< "edge:" << edgeI << " n0:" << n0 << endl;
+
+            bool same = true;
+
+            for (label i = 1; i < eFaces.size(); i++)
+            {
+                //const labelledTri& f = surf[i];
+                const Foam::vector& n = surf.faceNormals()[eFaces[i]];
+
+                //Pout<< "    mag(n&n0): " << mag(n&n0) << endl;
+
+                if (mag(n&n0) < minCos)
+                {
+                    same = false;
+                    break;
+                }
+            }
+
+            if (same)
+            {
+                edgeStat[edgeI] = surfaceFeatures::NONE;
+            }
+        }
+    }
+}
+
+
+Foam::surfaceFeatures::edgeStatus Foam::checkFlatRegionEdge
+(
+    const triSurface& surf,
+    const scalar tol,
+    const scalar includedAngle,
+    const label edgeI
+)
+{
+    const edge& e = surf.edges()[edgeI];
+    const labelList& eFaces = surf.edgeFaces()[edgeI];
+
+    // Bin according to normal
+
+    DynamicList<Foam::vector> normals(2);
+    DynamicList<labelList> bins(2);
+
+    forAll(eFaces, eFacei)
+    {
+        const Foam::vector& n = surf.faceNormals()[eFaces[eFacei]];
+
+        // Find the normal in normals
+        label index = -1;
+        forAll(normals, normalI)
+        {
+            if (mag(n&normals[normalI]) > (1-tol))
+            {
+                index = normalI;
+                break;
+            }
+        }
+
+        if (index != -1)
+        {
+            bins[index].append(eFacei);
+        }
+        else if (normals.size() >= 2)
+        {
+            // Would be third normal. Mark as feature.
+            //Pout<< "** at edge:" << surf.localPoints()[e[0]]
+            //    << surf.localPoints()[e[1]]
+            //    << " have normals:" << normals
+            //    << " and " << n << endl;
+            return surfaceFeatures::REGION;
+        }
+        else
+        {
+            normals.append(n);
+            bins.append(labelList(1, eFacei));
+        }
+    }
+
+
+    // Check resulting number of bins
+    if (bins.size() == 1)
+    {
+        // Note: should check here whether they are two sets of faces
+        // that are planar or indeed 4 faces al coming together at an edge.
+        //Pout<< "** at edge:"
+        //    << surf.localPoints()[e[0]]
+        //    << surf.localPoints()[e[1]]
+        //    << " have single normal:" << normals[0]
+        //    << endl;
+        return surfaceFeatures::NONE;
+    }
+    else
+    {
+        // Two bins. Check if normals make an angle
+
+        //Pout<< "** at edge:"
+        //    << surf.localPoints()[e[0]]
+        //    << surf.localPoints()[e[1]] << nl
+        //    << "    normals:" << normals << nl
+        //    << "    bins   :" << bins << nl
+        //    << endl;
+
+        if (includedAngle >= 0)
+        {
+            scalar minCos = Foam::cos(degToRad(180.0 - includedAngle));
+
+            forAll(eFaces, i)
+            {
+                const Foam::vector& ni = surf.faceNormals()[eFaces[i]];
+                for (label j=i+1; j<eFaces.size(); j++)
+                {
+                    const Foam::vector& nj = surf.faceNormals()[eFaces[j]];
+                    if (mag(ni & nj) < minCos)
+                    {
+                        //Pout<< "have sharp feature between normal:" << ni
+                        //    << " and " << nj << endl;
+
+                        // Is feature. Keep as region or convert to
+                        // feature angle? For now keep as region.
+                        return surfaceFeatures::REGION;
+                    }
+                }
+            }
+        }
+
+
+        // So now we have two normals bins but need to make sure both
+        // bins have the same regions in it.
+
+         // 1. store + or - region number depending
+        //    on orientation of triangle in bins[0]
+        const labelList& bin0 = bins[0];
+        labelList regionAndNormal(bin0.size());
+        forAll(bin0, i)
+        {
+            const labelledTri& t = surf.localFaces()[eFaces[bin0[i]]];
+            int dir = t.edgeDirection(e);
+
+            if (dir > 0)
+            {
+                regionAndNormal[i] = t.region()+1;
+            }
+            else if (dir == 0)
+            {
+                FatalErrorInFunction
+                    << exit(FatalError);
+            }
+            else
+            {
+                regionAndNormal[i] = -(t.region()+1);
+            }
+        }
+
+        // 2. check against bin1
+        const labelList& bin1 = bins[1];
+        labelList regionAndNormal1(bin1.size());
+        forAll(bin1, i)
+        {
+            const labelledTri& t = surf.localFaces()[eFaces[bin1[i]]];
+            int dir = t.edgeDirection(e);
+
+            label myRegionAndNormal;
+            if (dir > 0)
+            {
+                myRegionAndNormal = t.region()+1;
+            }
+            else
+            {
+                myRegionAndNormal = -(t.region()+1);
+            }
+
+            regionAndNormal1[i] = myRegionAndNormal;
+
+            label index = findIndex(regionAndNormal, -myRegionAndNormal);
+            if (index == -1)
+            {
+                // Not found.
+                //Pout<< "cannot find region " << myRegionAndNormal
+                //    << " in regions " << regionAndNormal << endl;
+
+                return surfaceFeatures::REGION;
+            }
+        }
+
+        // Passed all checks, two normal bins with the same contents.
+        //Pout<< "regionAndNormal:" << regionAndNormal << endl;
+        //Pout<< "myRegionAndNormal:" << regionAndNormal1 << endl;
+
+        return surfaceFeatures::NONE;
+    }
+}
+
+
+void Foam::extractCloseness
+(
+    const fileName &sFeatFileName,
+    const Time& runTime,
+    const triSurface &surf,
+    const bool writeVTK
+);
+
+
+void Foam::extractPointCloseness
+(
+    const fileName &sFeatFileName,
+    const Time& runTime,
+    const triSurface &surf,
+    const bool writeVTK
+);
+
+
+void Foam::writeStats(const extendedFeatureEdgeMesh& fem, Ostream& os)
+{
+    os  << "    points : " << fem.points().size() << nl
+        << "    of which" << nl
+        << "        convex             : "
+        << fem.concaveStart() << nl
+        << "        concave            : "
+        << (fem.mixedStart()-fem.concaveStart()) << nl
+        << "        mixed              : "
+        << (fem.nonFeatureStart()-fem.mixedStart()) << nl
+        << "        non-feature        : "
+        << (fem.points().size()-fem.nonFeatureStart()) << nl
+        << "    edges  : " << fem.edges().size() << nl
+        << "    of which" << nl
+        << "        external edges     : "
+        << fem.internalStart() << nl
+        << "        internal edges     : "
+        << (fem.flatStart()- fem.internalStart()) << nl
+        << "        flat edges         : "
+        << (fem.openStart()- fem.flatStart()) << nl
+        << "        open edges         : "
+        << (fem.multipleStart()- fem.openStart()) << nl
+        << "        multiply connected : "
+        << (fem.edges().size()- fem.multipleStart()) << endl;
+}
+
+
+// ************************************************************************* //

src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,9 +66,9 @@ bool Foam::triSurfaceSearch::checkUniqueHit
 
             if (pointFacei != currHit.index())
             {
-                forAll(hits, hI)
+                forAll(hits, hi)
                 {
-                    const pointIndexHit& hit = hits[hI];
+                    const pointIndexHit& hit = hits[hi];
 
                     if (hit.index() == pointFacei)
                     {
@@ -89,15 +89,15 @@ bool Foam::triSurfaceSearch::checkUniqueHit
 
         const labelList& edgeFaces = surface().edgeFaces()[edgeI];
 
-        forAll(edgeFaces, fI)
+        forAll(edgeFaces, fi)
         {
-            const label edgeFacei = edgeFaces[fI];
+            const label edgeFacei = edgeFaces[fi];
 
             if (edgeFacei != currHit.index())
             {
-                forAll(hits, hI)
+                forAll(hits, hi)
                 {
-                    const pointIndexHit& hit = hits[hI];
+                    const pointIndexHit& hit = hits[hi];
 
                     if (hit.index() == edgeFacei)
                     {
@@ -294,7 +294,7 @@ void Foam::triSurfaceSearch::findNearest
 
     forAll(samples, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findNearest
+        info[i] = octree.findNearest
         (
             samples[i],
             nearestDistSqr[i],
@@ -335,7 +335,7 @@ void Foam::triSurfaceSearch::findLine
 
     forAll(start, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findLine
+        info[i] = octree.findLine
         (
             start[i],
             end[i]
@@ -362,7 +362,7 @@ void Foam::triSurfaceSearch::findLineAny
 
     forAll(start, i)
     {
-        static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
+        info[i] = octree.findLineAny
         (
             start[i],
             end[i]
@@ -394,7 +394,7 @@ void Foam::triSurfaceSearch::findLineAll
 
     treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
 
-    forAll(start, pointi)
+    forAll(start, i)
     {
         hits.clear();
         shapeMask.clear();
@@ -404,25 +404,17 @@ void Foam::triSurfaceSearch::findLineAll
             // See if any intersection between pt and end
             pointIndexHit inter = octree.findLine
             (
-                start[pointi],
+                start[i],
-                end[pointi],
+                end[i],
                 allIntersectOp
             );
 
             if (inter.hit())
             {
-                vector lineVec = end[pointi] - start[pointi];
+                vector lineVec = end[i] - start[i];
                 lineVec /= mag(lineVec) + VSMALL;
 
-                if
+                if (checkUniqueHit(inter, hits, lineVec))
-                (
-                    checkUniqueHit
-                    (
-                        inter,
-                        hits,
-                        lineVec
-                    )
-                )
                 {
                     hits.append(inter);
                 }
@@ -435,7 +427,7 @@ void Foam::triSurfaceSearch::findLineAll
             }
         }
 
-        info[pointi].transfer(hits);
+        info[i].transfer(hits);
     }
 
     indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;

src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -170,7 +170,6 @@ public:
             const bool isNodeValues,
             const bool verbose = false
         ) const;
-
 };
 
 

src/triSurface/triSurface/triSurface.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1031,6 +1031,19 @@ Foam::triSurface Foam::triSurface::subsetMesh
 }
 
 
+Foam::faceList Foam::triSurface::faces() const
+{
+    faceList faces(size());
+
+    forAll(*this, facei)
+    {
+        faces[facei] = operator[](facei).triFaceFace();
+    }
+
+    return faces;
+}
+
+
 void Foam::triSurface::write
 (
     const fileName& name,

src/triSurface/triSurface/triSurface.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -393,6 +393,12 @@ public:
             ) const;
 
 
+        // Conversion
+
+            //- Return the list of triangles as a faceList
+            faceList faces() const;
+
+
         // Write
 
             //- Write to Ostream in simple FOAM format