ENH: surfaceCheck. New writeSets option. Fixes #717.

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,7 +142,7 @@ void writeZoning
               + '_'
               + surfFileNameBase
               + '.'
-              + vtkSurfaceWriter::typeName
+              + writer.type()
             )
         << " ..." << endl << endl;
 
@@ -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,8 +774,8 @@ int main(int argc, char *argv[])
                     }
                 }
 
-                nClose++;
+                if (nClose < outputThreshold)
-
+                {
                     if (edgei == -1)
                     {
                         Info<< "    close unconnected points "
@@ -633,6 +797,15 @@ int main(int argc, char *argv[])
                             << endl;
                     }
                 }
+                else if (nClose == outputThreshold)
+                {
+                    Info<< "Suppressing further warning messages."
+                        << " Use -outputThreshold to increase number"
+                        << " of warnings." << endl;
+                }
+
+                nClose++;
+            }
         }
 
         Info<< "Found " << nClose << " nearby points." << nl
@@ -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]);
-
+            openEdges.append(edgei);
-            nSingleEdges++;
         }
     }
 
-    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]);
             }
-
+            multipleEdges.append(edgei);
-            nMultEdges++;
         }
     }
     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 one face : " << openEdges.size() << endl
-            << "    connected to >2 faces : " << nMultEdges << 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");
+            const fileName openName
-
+            (
-            Info<< "Dumping conflicting face labels to " << str.name() << endl
+                surfFileName.lessExt()
-                << "Paste this into the input for surfaceSubset" << endl;
+              + "_open."
-
+              + edgeFormat
-            str << problemFaces;
+            );
+            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",

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;
@@ -2786,23 +2787,27 @@ bool Foam::triSurfaceTools::validTri
     forAll(f, fp)
     {
         if (f[fp] < 0 || f[fp] >= surf.points().size())
+        {
+            if (verbose)
             {
                 WarningInFunction
                     << "triangle " << facei << " vertices " << f
                     << " uses point indices outside point range 0.."
-                << surf.points().size()-1
+                    << surf.points().size()-1 << endl;
-                << endl;
+            }
             return false;
         }
     }
 
     if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
+    {
+        if (verbose)
         {
             WarningInFunction
                 << "triangle " << facei
                 << " uses non-unique vertices " << f
-            << " coords:" << f.points(surf.points())
+                << " coords:" << f.points(surf.points()) << endl;
-            << endl;
+        }
         return false;
     }
 
@@ -2831,13 +2836,15 @@ bool Foam::triSurfaceTools::validTri
          && (f[1] == nbrF[0] || f[1] == nbrF[1] || f[1] == nbrF[2])
          && (f[2] == nbrF[0] || f[2] == nbrF[1] || f[2] == nbrF[2])
         )
+        {
+            if (verbose)
             {
                 WarningInFunction
                     << "triangle " << facei << " vertices " << f
                     << " has the same vertices as triangle " << nbrFacei
                     << " vertices " << nbrF
-                << " coords:" << f.points(surf.points())
+                    << " coords:" << f.points(surf.points()) << endl;
-                << endl;
+            }
 
             return false;
         }
@@ -2850,19 +2857,22 @@ bool Foam::triSurfaceTools::validTri
 bool Foam::triSurfaceTools::validTri
 (
     const MeshedSurface<face>& surf,
-    const label facei
+    const label facei,
+    const bool verbose
 )
 {
     typedef face FaceType;
     const FaceType& f = surf[facei];
 
     if (f.size() != 3)
+    {
+        if (verbose)
         {
             WarningInFunction
                 << "face " << facei
                 << " is not a triangle, it has " << f.size()
-            << " indices"
+                << " indices" << endl;
-            << endl;
+        }
         return false;
     }
 
@@ -2870,23 +2880,27 @@ bool Foam::triSurfaceTools::validTri
     forAll(f, fp)
     {
         if (f[fp] < 0 || f[fp] >= surf.points().size())
+        {
+            if (verbose)
             {
                 WarningInFunction
                     << "triangle " << facei << " vertices " << f
                     << " uses point indices outside point range 0.."
-                << surf.points().size()-1
+                    << surf.points().size()-1 << endl;
-                << endl;
+            }
             return false;
         }
     }
 
     if (f[0] == f[1] || f[0] == f[2] || f[1] == f[2])
+    {
+        if (verbose)
         {
             WarningInFunction
                 << "triangle " << facei
                 << " uses non-unique vertices " << f
-            << " coords:" << f.points(surf.points())
+                << " coords:" << f.points(surf.points()) << endl;
-            << endl;
+        }
         return false;
     }
 
@@ -2915,14 +2929,15 @@ bool Foam::triSurfaceTools::validTri
          && (f[1] == nbrF[0] || f[1] == nbrF[1] || f[1] == nbrF[2])
          && (f[2] == nbrF[0] || f[2] == nbrF[1] || f[2] == nbrF[2])
         )
+        {
+            if (verbose)
             {
                 WarningInFunction
                     << "triangle " << facei << " vertices " << f
                     << " has the same vertices as triangle " << nbrFacei
                     << " vertices " << nbrF
-                << " coords:" << f.points(surf.points())
+                    << " coords:" << f.points(surf.points()) << endl;
-                << endl;
+            }
-
             return false;
         }
     }

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