surfaceFeatureExtract: Further refactoring to move the extract surface closeness functionality

to triSurfaceMesh

src/meshTools/Make/files

@@ -114,6 +114,7 @@ searchableSurfaces/searchableSurfaces/searchableSurfaces.C
 searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C
 searchableSurfaces/searchableSurfaceWithGaps/searchableSurfaceWithGaps.C
 searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
+searchableSurfaces/triSurfaceMesh/extractCloseness.C
 searchableSurfaces/closedTriSurfaceMesh/closedTriSurfaceMesh.C
 searchableSurfaces/searchableExtrudedCircle/searchableExtrudedCircle.C
 

src/meshTools/searchableSurfaces/triSurfaceMesh/extractCloseness.C

@@ -0,0 +1,425 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "triSurfaceMesh.H"
+#include "triSurfaceFields.H"
+#include "meshTools.H"
+#include "Time.H"
+#include "unitConversion.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::triSurfaceMesh::drawHitProblem
+(
+    const label fi,
+    const point& start,
+    const point& p,
+    const point& end,
+    const pointIndexHitList& hitInfo
+) const
+{
+    const List<labelledTri>& tris = *this;
+    const pointField& points = this->points();
+
+    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, points[tris[fi][0]]);
+    meshTools::writeOBJ(Info, points[tris[fi][1]]);
+    meshTools::writeOBJ(Info, points[tris[fi][2]]);
+
+    Info<< "f 4 5 6" << endl;
+
+    forAll(hitInfo, hi)
+    {
+        label hfi = hitInfo[hi].index();
+
+        meshTools::writeOBJ(Info, points[tris[hfi][0]]);
+        meshTools::writeOBJ(Info, points[tris[hfi][1]]);
+        meshTools::writeOBJ(Info, points[tris[hfi][2]]);
+
+        Info<< "f "
+            << 3*hi + 7 << " "
+            << 3*hi + 8 << " "
+            << 3*hi + 9
+            << endl;
+    }
+}
+
+
+void Foam::triSurfaceMesh::processHit
+(
+    scalar& internalCloseness,
+    scalar& externalCloseness,
+    const scalar internalToleranceCosAngle,
+    const scalar externalToleranceCosAngle,
+    const label fi,
+    const point& start,
+    const point& p,
+    const point& end,
+    const vector& normal,
+    const vectorField& normals,
+    const pointIndexHitList& hitInfo
+) const
+{
+    if (hitInfo.size() < 1)
+    {
+        drawHitProblem(fi, start, p, end, hitInfo);
+    }
+    else if (hitInfo.size() == 1)
+    {
+        if (!hitInfo[0].hit())
+        {
+        }
+        else if (hitInfo[0].index() != fi)
+        {
+            drawHitProblem(fi, 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, 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())
+                );
+            }
+        }
+    }
+}
+
+
+Foam::Pair<Foam::tmp<Foam::triSurfaceScalarField>>
+Foam::triSurfaceMesh::extractCloseness
+(
+    const scalar internalAngleTolerance,
+    const scalar externalAngleTolerance
+) const
+{
+    const scalar internalToleranceCosAngle
+    (
+        cos(degToRad(180 - internalAngleTolerance))
+    );
+
+    const scalar externalToleranceCosAngle
+    (
+        cos(degToRad(180 - externalAngleTolerance))
+    );
+
+    const Time& runTime = objectRegistry::time();
+
+    // Prepare start and end points for intersection tests
+
+    const vectorField& normals = faceNormals();
+
+    const scalar span = bounds().mag();
+
+    const pointField start(faceCentres() - span*normals);
+    const pointField end(faceCentres() + span*normals);
+    const pointField& faceCentres = this->faceCentres();
+
+    List<List<pointIndexHit>> allHitinfo;
+
+    // Find all intersections (in order)
+    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],
+            internalToleranceCosAngle,
+            externalToleranceCosAngle,
+            fi,
+            start[fi],
+            faceCentres[fi],
+            end[fi],
+            normals[fi],
+            normals,
+            hitInfo
+        );
+    }
+
+    return Pair<tmp<triSurfaceScalarField>>
+    (
+        tmp<triSurfaceScalarField>
+        (
+            new triSurfaceScalarField
+            (
+                IOobject
+                (
+                    objectRegistry::name() + ".internalCloseness",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime
+                ),
+                *this,
+                dimLength,
+                internalCloseness
+            )
+        ),
+
+        tmp<triSurfaceScalarField>
+        (
+            new triSurfaceScalarField
+            (
+                IOobject
+                (
+                    objectRegistry::name() + ".externalCloseness",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime
+                ),
+                *this,
+                dimLength,
+                externalCloseness
+            )
+        )
+    );
+}
+
+
+Foam::Pair<Foam::tmp<Foam::triSurfacePointScalarField>>
+Foam::triSurfaceMesh::extractPointCloseness
+(
+    const scalar internalAngleTolerance,
+    const scalar externalAngleTolerance
+) const
+{
+    const scalar internalToleranceCosAngle
+    (
+        cos(degToRad(180 - internalAngleTolerance))
+    );
+
+    const scalar externalToleranceCosAngle
+    (
+        cos(degToRad(180 - externalAngleTolerance))
+    );
+
+    const Time& runTime = objectRegistry::time();
+
+    // Prepare start and end points for intersection tests
+
+    const pointField& points = this->points();
+    const labelList& meshPoints = this->meshPoints();
+    const pointField& faceCentres = this->faceCentres();
+    const vectorField& normals = this->faceNormals();
+    const labelListList& pointFaces = this->pointFaces();
+
+    const scalar span = 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<pointIndexHitList> allHitinfo;
+
+    // Find all intersections (in order)
+    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 pointIndexHitList& hitInfo = allHitinfo[i];
+
+            processHit
+            (
+                internalCloseness[pi],
+                externalCloseness[pi],
+                internalToleranceCosAngle,
+                externalToleranceCosAngle,
+                fi,
+                start[i],
+                facePoints[i],
+                end[i],
+                normals[fi],
+                normals,
+                hitInfo
+            );
+
+            i++;
+        }
+    }
+
+    return Pair<tmp<triSurfacePointScalarField>>
+    (
+        tmp<triSurfacePointScalarField>
+        (
+            new triSurfacePointScalarField
+            (
+                IOobject
+                (
+                    objectRegistry::name() + ".internalPointCloseness",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime
+                ),
+                *this,
+                dimLength,
+                internalCloseness
+            )
+        ),
+
+        tmp<triSurfacePointScalarField>
+        (
+            new triSurfacePointScalarField
+            (
+                IOobject
+                (
+                    objectRegistry::name() + ".externalPointCloseness",
+                    runTime.constant(),
+                    "triSurface",
+                    runTime
+                ),
+                *this,
+                dimLength,
+                externalCloseness
+            )
+        )
+    );
+}
+
+
+// ************************************************************************* //

src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C

@@ -794,6 +794,37 @@ void Foam::triSurfaceMesh::getNormal
 }
 
 
+void Foam::triSurfaceMesh::getVolumeType
+(
+    const pointField& points,
+    List<volumeType>& volType
+) const
+{
+    volType.setSize(points.size());
+
+    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
+    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
+
+    forAll(points, pointi)
+    {
+        const point& pt = points[pointi];
+
+        if (!tree().bb().contains(pt))
+        {
+            // Have to calculate directly as outside the octree
+            volType[pointi] = tree().shapes().getVolumeType(tree(), pt);
+        }
+        else
+        {
+            // - use cached volume type per each tree node
+            volType[pointi] = tree().getVolumeType(pt);
+        }
+    }
+
+    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
+}
+
+
 void Foam::triSurfaceMesh::setField(const labelList& values)
 {
     autoPtr<triSurfaceLabelField> fldPtr
@@ -846,37 +877,6 @@ void Foam::triSurfaceMesh::getField
 }
 
 
-void Foam::triSurfaceMesh::getVolumeType
-(
-    const pointField& points,
-    List<volumeType>& volType
-) const
-{
-    volType.setSize(points.size());
-
-    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
-    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
-
-    forAll(points, pointi)
-    {
-        const point& pt = points[pointi];
-
-        if (!tree().bb().contains(pt))
-        {
-            // Have to calculate directly as outside the octree
-            volType[pointi] = tree().shapes().getVolumeType(tree(), pt);
-        }
-        else
-        {
-            // - use cached volume type per each tree node
-            volType[pointi] = tree().getVolumeType(pt);
-        }
-    }
-
-    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
-}
-
-
 bool Foam::triSurfaceMesh::writeObject
 (
     IOstream::streamFormat fmt,

src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H

@@ -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
@@ -51,6 +51,8 @@ SourceFiles
 #include "EdgeMap.H"
 #include "triSurface.H"
 #include "triSurfaceRegionSearch.H"
+#include "triSurfaceFieldsFwd.H"
+#include "pointIndexHitList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -131,6 +133,30 @@ class triSurfaceMesh
             DynamicList<pointIndexHit, 1, 1>& hits
         );
 
+        void drawHitProblem
+        (
+            const label fi,
+            const point& start,
+            const point& p,
+            const point& end,
+            const pointIndexHitList& hitInfo
+        ) const;
+
+        void processHit
+        (
+            scalar& internalCloseness,
+            scalar& externalCloseness,
+            const scalar internalToleranceCosAngle,
+            const scalar externalToleranceCosAngle,
+            const label fi,
+            const point& start,
+            const point& p,
+            const point& end,
+            const vector& normal,
+            const vectorField& normals,
+            const pointIndexHitList& hitInfo
+        ) const;
+
         //- Disallow default bitwise copy construct
         triSurfaceMesh(const triSurfaceMesh&);
 
@@ -290,6 +316,18 @@ public:
             //  indices) get the specified field. Misses do not get set.
             virtual void getField(const List<pointIndexHit>&, labelList&) const;
 
+            Pair<tmp<triSurfaceScalarField>> extractCloseness
+            (
+                const scalar internalAngleTolerance = 80,
+                const scalar externalAngleTolerance = 10
+            ) const;
+
+            Pair<tmp<triSurfacePointScalarField>> extractPointCloseness
+            (
+                const scalar internalAngleTolerance = 80,
+                const scalar externalAngleTolerance = 10
+            ) const;
+
 
         // regIOobject implementation
 

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C

@@ -1499,4 +1499,167 @@ void Foam::surfaceFeatures::operator=(const surfaceFeatures& rhs)
 }
 
 
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+Foam::surfaceFeatures::edgeStatus Foam::checkNonManifoldEdge
+(
+    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;
+            }
+        }
+
+        return surfaceFeatures::NONE;
+    }
+}
+
+
 // ************************************************************************* //

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.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
@@ -32,13 +32,11 @@ Description
         externalStart_ .. internalStart_-1  : external edges
         internalStart_ .. size-1            : internal edges
 
-
     NOTE: angle is included angle, not feature angle and is in degrees.
     The included angle is the smallest angle between two planes. For coplanar
     faces it is 180, for straight angles it is 90. To pick up straight edges
     only use included angle of 91 degrees
 
-
 SourceFiles
     surfaceFeatures.C
 
@@ -418,11 +416,22 @@ public:
     // Member Operators
 
         void operator=(const surfaceFeatures&);
-
-
 };
 
 
+//- 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 checkNonManifoldEdge
+(
+    const triSurface& surf,
+    const scalar tol,
+    const scalar includedAngle,
+    const label edgei
+);
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam