From 9ee6036bb29455b6e5044ce8b3bda8a03a8d3b84 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Thu, 12 Jul 2018 13:33:28 +0100
Subject: [PATCH] ENH: surfaceCheck. New writeSets option. Fixes #717.

---
 .../surface/surfaceCheck/surfaceCheck.C       | 286 +++++++++++++++---
 .../triSurfaceTools/triSurfaceTools.C         |  95 +++---
 .../triSurfaceTools/triSurfaceTools.H         |  14 +-
 3 files changed, 311 insertions(+), 84 deletions(-)

diff --git a/applications/utilities/surface/surfaceCheck/surfaceCheck.C b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
index 1b76d0a8add..0b3acc607ec 100644
--- a/applications/utilities/surface/surfaceCheck/surfaceCheck.C
+++ b/applications/utilities/surface/surfaceCheck/surfaceCheck.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2016-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,6 +64,9 @@ Usage
 #include "SortableList.H"
 #include "PatchTools.H"
 #include "vtkSurfaceWriter.H"
+#include "functionObject.H"
+#include "DynamicField.H"
+#include "edgeMesh.H"
 
 using namespace Foam;
 
@@ -123,6 +126,7 @@ labelList countBins
 
 void writeZoning
 (
+    const surfaceWriter& writer,
     const triSurface& surf,
     const labelList& faceZone,
     const word& fieldName,
@@ -138,9 +142,9 @@ void writeZoning
               + '_'
               + surfFileNameBase
               + '.'
-              + vtkSurfaceWriter::typeName
+              + writer.type()
             )
-        << "..." << endl << endl;
+        << " ..." << endl << endl;
 
     // Convert data
     scalarField scalarFaceZone(faceZone.size());
@@ -154,7 +158,7 @@ void writeZoning
         faces[i] = surf[i];
     }
 
-    vtkSurfaceWriter().write
+    writer.write
     (
         surfFilePath,
         surfFileNameBase,
@@ -275,6 +279,40 @@ void syncEdges(const triSurface& p, boolList& isMarkedEdge)
 }
 
 
+void writeEdgeSet
+(
+    const word& setName,
+    const triSurface& surf,
+    const labelUList& edgeSet
+)
+{
+    // Get compact edge mesh
+    labelList pointToCompact(surf.nPoints(), -1);
+    DynamicField<point> compactPoints(edgeSet.size());
+    DynamicList<edge> compactEdges(edgeSet.size());
+    for (label edgei : edgeSet)
+    {
+        const edge& e = surf.edges()[edgei];
+        edge compactEdge(-1, -1);
+        forAll(e, ep)
+        {
+            label& compacti = pointToCompact[e[ep]];
+            if (compacti == -1)
+            {
+                compacti = compactPoints.size();
+                label pointi = surf.meshPoints()[e[ep]];
+                compactPoints.append(surf.points()[pointi]);
+            }
+            compactEdge[ep] = compacti;
+        }
+        compactEdges.append(compactEdge);
+    }
+
+    edgeMesh eMesh(std::move(compactPoints), std::move(compactEdges));
+    eMesh.write(setName);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -309,6 +347,13 @@ int main(int argc, char *argv[])
         "upper limit on the number of files written."
         " Default is 10, using 0 suppresses file writing."
     );
+    argList::addOption
+    (
+        "writeSets",
+        "surfaceFormat",
+        "reconstruct and write problem triangles/edges in selected format"
+    );
+
 
     argList args(argc, argv);
 
@@ -317,6 +362,25 @@ int main(int argc, char *argv[])
     const bool splitNonManifold = args.found("splitNonManifold");
     const label outputThreshold =
         args.lookupOrDefault("outputThreshold", 10);
+    word surfaceFormat;
+    const bool writeSets = args.optionReadIfPresent("writeSets", surfaceFormat);
+
+    autoPtr<surfaceWriter> surfWriter;
+    word edgeFormat;
+    if (writeSets)
+    {
+        surfWriter = surfaceWriter::New(surfaceFormat);
+        // Option1: hard-coded format
+        edgeFormat = "obj";
+        //// Option2: same type as surface format. Problem is e.g. .obj format
+        ////          is not a surfaceWriter format.
+        //edgeFormat = surfaceFormat;
+        //if (!edgeMesh::canWriteType(edgeFormat))
+        //{
+        //    edgeFormat = "obj";
+        //}
+    }
+
 
     Info<< "Reading surface from " << surfFileName << " ..." << nl << endl;
 
@@ -414,7 +478,8 @@ int main(int argc, char *argv[])
 
         forAll(surf, facei)
         {
-            if (!triSurfaceTools::validTri(surf, facei))
+            // Check silently
+            if (!triSurfaceTools::validTri(surf, facei, false))
             {
                 illegalFaces.append(facei);
             }
@@ -425,8 +490,71 @@ int main(int argc, char *argv[])
             Info<< "Surface has " << illegalFaces.size()
                 << " illegal triangles." << endl;
 
-            if (outputThreshold > 0)
+            if (surfWriter.valid())
+            {
+                boolList isIllegalFace(surf.size(), false);
+                UIndirectList<bool>(isIllegalFace, illegalFaces) = true;
+
+                labelList pointMap;
+                labelList faceMap;
+                triSurface subSurf
+                (
+                    surf.subsetMesh
+                    (
+                        isIllegalFace,
+                        pointMap,
+                        faceMap
+                    )
+                );
+
+                const fileName qualityName
+                (
+                    surfFilePath
+                  / "illegal"
+                  + '_'
+                  + surfFileNameBase
+                  + '.'
+                  + surfWriter().type()
+                );
+                Info<< "Writing illegal triangles to "
+                    << qualityName << " ..." << endl << endl;
+
+                // Convert data
+                faceList faces(subSurf.size());
+                forAll(subSurf, i)
+                {
+                    faces[i] = subSurf[i];
+                }
+
+                surfWriter().write
+                (
+                    surfFilePath,
+                    surfFileNameBase,
+                    meshedSurfRef
+                    (
+                        subSurf.points(),
+                        faces
+                    ),
+                    "illegal",
+                    scalarField(subSurf.size(), 0.0),
+                    false               // face based data
+                );
+            }
+            else if (outputThreshold > 0)
             {
+                forAll(illegalFaces, i)
+                {
+                    // Force warning message
+                    triSurfaceTools::validTri(surf, illegalFaces[i], true);
+                    if (i >= outputThreshold)
+                    {
+                        Info<< "Suppressing further warning messages."
+                            << " Use -outputThreshold to increase number"
+                            << " of warnings." << endl;
+                        break;
+                    }
+                }
+
                 OFstream str("illegalFaces");
                 Info<< "Dumping conflicting face labels to " << str.name()
                     << endl
@@ -502,12 +630,48 @@ int main(int argc, char *argv[])
         if (triQ[minIndex] < SMALL)
         {
             WarningInFunction
+                << "Minimum triangle quality is "
                 << triQ[minIndex] << ". This might give problems in"
                 << " self-intersection testing later on." << endl;
         }
 
         // Dump for subsetting
-        if (outputThreshold > 0)
+        if (surfWriter.valid())
+        {
+            const fileName qualityName
+            (
+                surfFilePath
+              / "quality"
+              + '_'
+              + surfFileNameBase
+              + '.'
+              + surfWriter().type()
+            );
+            Info<< "Writing triangle-quality to "
+                << qualityName << " ..." << endl << endl;
+
+            // Convert data
+            faceList faces(surf.size());
+            forAll(surf, i)
+            {
+                faces[i] = surf[i];
+            }
+
+            surfWriter().write
+            (
+                surfFilePath,
+                surfFileNameBase,
+                meshedSurfRef
+                (
+                    surf.points(),
+                    faces
+                ),
+                "quality",
+                triQ,
+                false               // face based data
+            );
+        }
+        else if (outputThreshold > 0)
         {
             DynamicList<label> problemFaces(surf.size()/100+1);
 
@@ -610,28 +774,37 @@ int main(int argc, char *argv[])
                     }
                 }
 
-                nClose++;
-
-                if (edgei == -1)
+                if (nClose < outputThreshold)
                 {
-                    Info<< "    close unconnected points "
-                        << pti << ' ' << localPoints[pti]
-                        << " and " << prevPti << ' '
-                        << localPoints[prevPti]
-                        << " distance:"
-                        << mag(localPoints[pti] - localPoints[prevPti])
-                        << endl;
+                    if (edgei == -1)
+                    {
+                        Info<< "    close unconnected points "
+                            << pti << ' ' << localPoints[pti]
+                            << " and " << prevPti << ' '
+                            << localPoints[prevPti]
+                            << " distance:"
+                            << mag(localPoints[pti] - localPoints[prevPti])
+                            << endl;
+                    }
+                    else
+                    {
+                        Info<< "    small edge between points "
+                            << pti << ' ' << localPoints[pti]
+                            << " and " << prevPti << ' '
+                            << localPoints[prevPti]
+                            << " distance:"
+                            << mag(localPoints[pti] - localPoints[prevPti])
+                            << endl;
+                    }
                 }
-                else
+                else if (nClose == outputThreshold)
                 {
-                    Info<< "    small edge between points "
-                        << pti << ' ' << localPoints[pti]
-                        << " and " << prevPti << ' '
-                        << localPoints[prevPti]
-                        << " distance:"
-                        << mag(localPoints[pti] - localPoints[prevPti])
-                        << endl;
+                    Info<< "Suppressing further warning messages."
+                        << " Use -outputThreshold to increase number"
+                        << " of warnings." << endl;
                 }
+
+                nClose++;
             }
         }
 
@@ -645,10 +818,11 @@ int main(int argc, char *argv[])
     // ~~~~~~~~~~~~~~
 
     DynamicList<label> problemFaces(surf.size()/100 + 1);
+    DynamicList<label> openEdges(surf.nEdges()/100 + 1);
+    DynamicList<label> multipleEdges(surf.nEdges()/100 + 1);
 
     const labelListList& edgeFaces = surf.edgeFaces();
 
-    label nSingleEdges = 0;
     forAll(edgeFaces, edgei)
     {
         const labelList& myFaces = edgeFaces[edgei];
@@ -656,12 +830,10 @@ int main(int argc, char *argv[])
         if (myFaces.size() == 1)
         {
             problemFaces.append(myFaces[0]);
-
-            nSingleEdges++;
+            openEdges.append(edgei);
         }
     }
 
-    label nMultEdges = 0;
     forAll(edgeFaces, edgei)
     {
         const labelList& myFaces = edgeFaces[edgei];
@@ -672,30 +844,43 @@ int main(int argc, char *argv[])
             {
                 problemFaces.append(myFaces[myFacei]);
             }
-
-            nMultEdges++;
+            multipleEdges.append(edgei);
         }
     }
     problemFaces.shrink();
 
-    if ((nSingleEdges != 0) || (nMultEdges != 0))
+    if (openEdges.size() || multipleEdges.size())
     {
         Info<< "Surface is not closed since not all edges ("
             << edgeFaces.size() << ") connected to "
             << "two faces:" << endl
-            << "    connected to one face : " << nSingleEdges << endl
-            << "    connected to >2 faces : " << nMultEdges << endl;
+            << "    connected to one face : " << openEdges.size() << endl
+            << "    connected to >2 faces : " << multipleEdges.size() << endl;
 
         Info<< "Conflicting face labels:" << problemFaces.size() << endl;
 
-        if (outputThreshold > 0)
+        if (!edgeFormat.empty() && openEdges.size())
         {
-            OFstream str("problemFaces");
-
-            Info<< "Dumping conflicting face labels to " << str.name() << endl
-                << "Paste this into the input for surfaceSubset" << endl;
-
-            str << problemFaces;
+            const fileName openName
+            (
+                surfFileName.lessExt()
+              + "_open."
+              + edgeFormat
+            );
+            Info<< "Writing open edges to " << openName << " ..." << endl;
+            writeEdgeSet(openName, surf, openEdges);
+        }
+        if (!edgeFormat.empty() && multipleEdges.size())
+        {
+            const fileName multName
+            (
+                surfFileName.lessExt()
+              + "_multiply."
+              + edgeFormat
+            );
+            Info<< "Writing multiply connected edges to "
+                << multName << " ..." << endl;
+            writeEdgeSet(multName, surf, multipleEdges);
         }
     }
     else
@@ -732,7 +917,19 @@ int main(int argc, char *argv[])
         {
             Info<< "Splitting surface into parts ..." << endl << endl;
 
-            writeZoning(surf, faceZone, "zone", surfFilePath, surfFileNameBase);
+            if (!surfWriter.valid())
+            {
+                surfWriter.reset(new vtkSurfaceWriter());
+            }
+            writeZoning
+            (
+                surfWriter(),
+                surf,
+                faceZone,
+                "zone",
+                surfFilePath,
+                surfFileNameBase
+            );
 
             if (numZones > outputThreshold)
             {
@@ -785,8 +982,13 @@ int main(int argc, char *argv[])
 
         if (outputThreshold > 0)
         {
+            if (!surfWriter.valid())
+            {
+                surfWriter.reset(new vtkSurfaceWriter());
+            }
             writeZoning
             (
+                surfWriter(),
                 surf,
                 normalZone,
                 "normal",
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 7a47b6393a9..5b656338921 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -2776,7 +2776,8 @@ void Foam::triSurfaceTools::calcInterpolationWeights
 bool Foam::triSurfaceTools::validTri
 (
     const triSurface& surf,
-    const label facei
+    const label facei,
+    const bool verbose
 )
 {
     typedef labelledTri FaceType;
@@ -2787,22 +2788,26 @@ bool Foam::triSurfaceTools::validTri
     {
         if (f[fp] < 0 || f[fp] >= surf.points().size())
         {
-            WarningInFunction
-                << "triangle " << facei << " vertices " << f
-                << " uses point indices outside point range 0.."
-                << surf.points().size()-1
-                << endl;
+            if (verbose)
+            {
+                WarningInFunction
+                    << "triangle " << facei << " vertices " << f
+                    << " uses point indices outside point range 0.."
+                    << surf.points().size()-1 << endl;
+            }
             return false;
         }
     }
 
     if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
     {
-        WarningInFunction
-            << "triangle " << facei
-            << " uses non-unique vertices " << f
-            << " coords:" << f.points(surf.points())
-            << endl;
+        if (verbose)
+        {
+            WarningInFunction
+                << "triangle " << facei
+                << " uses non-unique vertices " << f
+                << " coords:" << f.points(surf.points()) << endl;
+        }
         return false;
     }
 
@@ -2832,12 +2837,14 @@ bool Foam::triSurfaceTools::validTri
          && (f[2] == nbrF[0] || f[2] == nbrF[1] || f[2] == nbrF[2])
         )
         {
-            WarningInFunction
-                << "triangle " << facei << " vertices " << f
-                << " has the same vertices as triangle " << nbrFacei
-                << " vertices " << nbrF
-                << " coords:" << f.points(surf.points())
-                << endl;
+            if (verbose)
+            {
+                WarningInFunction
+                    << "triangle " << facei << " vertices " << f
+                    << " has the same vertices as triangle " << nbrFacei
+                    << " vertices " << nbrF
+                    << " coords:" << f.points(surf.points()) << endl;
+            }
 
             return false;
         }
@@ -2850,7 +2857,8 @@ bool Foam::triSurfaceTools::validTri
 bool Foam::triSurfaceTools::validTri
 (
     const MeshedSurface<face>& surf,
-    const label facei
+    const label facei,
+    const bool verbose
 )
 {
     typedef face FaceType;
@@ -2858,11 +2866,13 @@ bool Foam::triSurfaceTools::validTri
 
     if (f.size() != 3)
     {
-        WarningInFunction
-            << "face " << facei
-            << " is not a triangle, it has " << f.size()
-            << " indices"
-            << endl;
+        if (verbose)
+        {
+            WarningInFunction
+                << "face " << facei
+                << " is not a triangle, it has " << f.size()
+                << " indices" << endl;
+        }
         return false;
     }
 
@@ -2871,22 +2881,26 @@ bool Foam::triSurfaceTools::validTri
     {
         if (f[fp] < 0 || f[fp] >= surf.points().size())
         {
-            WarningInFunction
-                << "triangle " << facei << " vertices " << f
-                << " uses point indices outside point range 0.."
-                << surf.points().size()-1
-                << endl;
+            if (verbose)
+            {
+                WarningInFunction
+                    << "triangle " << facei << " vertices " << f
+                    << " uses point indices outside point range 0.."
+                    << surf.points().size()-1 << endl;
+            }
             return false;
         }
     }
 
     if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
     {
-        WarningInFunction
-            << "triangle " << facei
-            << " uses non-unique vertices " << f
-            << " coords:" << f.points(surf.points())
-            << endl;
+        if (verbose)
+        {
+            WarningInFunction
+                << "triangle " << facei
+                << " uses non-unique vertices " << f
+                << " coords:" << f.points(surf.points()) << endl;
+        }
         return false;
     }
 
@@ -2916,13 +2930,14 @@ bool Foam::triSurfaceTools::validTri
          && (f[2] == nbrF[0] || f[2] == nbrF[1] || f[2] == nbrF[2])
         )
         {
-            WarningInFunction
-                << "triangle " << facei << " vertices " << f
-                << " has the same vertices as triangle " << nbrFacei
-                << " vertices " << nbrF
-                << " coords:" << f.points(surf.points())
-                << endl;
-
+            if (verbose)
+            {
+                WarningInFunction
+                    << "triangle " << facei << " vertices " << f
+                    << " has the same vertices as triangle " << nbrFacei
+                    << " vertices " << nbrF
+                    << " coords:" << f.points(surf.points()) << endl;
+            }
             return false;
         }
     }
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
index a761186b903..4ca5de6c963 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.H
@@ -617,10 +617,20 @@ public:
     // Surface checking functionality
 
         //- Check single triangle for (topological) validity
-        static bool validTri(const triSurface&, const label facei);
+        static bool validTri
+        (
+            const triSurface&,
+            const label facei,
+            const bool verbose = true
+        );
 
         //- Check single triangle for (topological) validity
-        static bool validTri(const MeshedSurface<face>&, const label facei);
+        static bool validTri
+        (
+            const MeshedSurface<face>&,
+            const label facei,
+            const bool verbose = true
+        );
 
 
     // Tracking
-- 
GitLab