checkMesh: Added option to write sets

- the checking for point-connected multiple-regions now also writes the
    conflicting points to a pointSet
  - with the -writeSets option it now also reconstructs & writes pointSets

applications/utilities/mesh/manipulation/checkMesh/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/surfMesh/lnInclude \

applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C

@@ -8,9 +8,11 @@
 #include "wedgePolyPatch.H"
 #include "unitConversion.H"
 #include "polyMeshTetDecomposition.H"
-#include "surfaceWriter.H"
+
-#include "checkTools.H"
 #include "vtkSurfaceWriter.H"
+#include "writer.H"
+
+#include "checkTools.H"
 #include "cyclicAMIPolyPatch.H"
 #include "Time.H"
 
@@ -481,7 +483,8 @@ Foam::label Foam::checkGeometry
 (
     const polyMesh& mesh,
     const bool allGeometry,
-    const autoPtr<surfaceWriter>& writer
+    const autoPtr<surfaceWriter>& surfWriter,
+    const autoPtr<writer<scalar>>& setWriter
 )
 {
     label noFailedChecks = 0;
@@ -538,6 +541,10 @@ Foam::label Foam::checkGeometry
                     << nonAlignedPoints.name() << endl;
                 nonAlignedPoints.instance() = mesh.pointsInstance();
                 nonAlignedPoints.write();
+                if (setWriter.valid())
+                {
+                    mergeAndWrite(setWriter, nonAlignedPoints);
+                }
             }
         }
     }
@@ -568,9 +575,9 @@ Foam::label Foam::checkGeometry
                     << " non closed cells to set " << cells.name() << endl;
                 cells.instance() = mesh.pointsInstance();
                 cells.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), cells);
+                    mergeAndWrite(surfWriter(), cells);
                 }
             }
         }
@@ -584,9 +591,9 @@ Foam::label Foam::checkGeometry
                 << aspectCells.name() << endl;
             aspectCells.instance() = mesh.pointsInstance();
             aspectCells.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), aspectCells);
+                mergeAndWrite(surfWriter(), aspectCells);
             }
         }
     }
@@ -605,9 +612,9 @@ Foam::label Foam::checkGeometry
                     << " zero area faces to set " << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), faces);
+                    mergeAndWrite(surfWriter(), faces);
                 }
             }
         }
@@ -627,9 +634,9 @@ Foam::label Foam::checkGeometry
                     << " zero volume cells to set " << cells.name() << endl;
                 cells.instance() = mesh.pointsInstance();
                 cells.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), cells);
+                    mergeAndWrite(surfWriter(), cells);
                 }
             }
         }
@@ -650,9 +657,9 @@ Foam::label Foam::checkGeometry
                 << " non-orthogonal faces to set " << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), faces);
+                mergeAndWrite(surfWriter(), faces);
             }
         }
     }
@@ -672,9 +679,9 @@ Foam::label Foam::checkGeometry
                     << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), faces);
+                    mergeAndWrite(surfWriter(), faces);
                 }
             }
         }
@@ -694,9 +701,9 @@ Foam::label Foam::checkGeometry
                     << " skew faces to set " << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), faces);
+                    mergeAndWrite(surfWriter(), faces);
                 }
             }
         }
@@ -718,9 +725,9 @@ Foam::label Foam::checkGeometry
                     << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), faces);
+                    mergeAndWrite(surfWriter(), faces);
                 }
             }
         }
@@ -751,9 +758,9 @@ Foam::label Foam::checkGeometry
                     << "decomposition tets to set " << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), faces);
+                    mergeAndWrite(surfWriter(), faces);
                 }
             }
         }
@@ -776,6 +783,10 @@ Foam::label Foam::checkGeometry
                     << endl;
                 points.instance() = mesh.pointsInstance();
                 points.write();
+                if (setWriter.valid())
+                {
+                    mergeAndWrite(setWriter, points);
+                }
             }
         }
 
@@ -795,6 +806,10 @@ Foam::label Foam::checkGeometry
                     << " apart) points to set " << nearPoints.name() << endl;
                 nearPoints.instance() = mesh.pointsInstance();
                 nearPoints.write();
+                if (setWriter.valid())
+                {
+                    mergeAndWrite(setWriter, nearPoints);
+                }
             }
         }
     }
@@ -815,9 +830,9 @@ Foam::label Foam::checkGeometry
                     << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), faces);
+                    mergeAndWrite(surfWriter(), faces);
                 }
             }
         }
@@ -838,9 +853,9 @@ Foam::label Foam::checkGeometry
                     << " warped faces to set " << faces.name() << endl;
                 faces.instance() = mesh.pointsInstance();
                 faces.write();
-                if (writer.valid())
+                if (surfWriter.valid())
                 {
-                    mergeAndWrite(writer(), faces);
+                    mergeAndWrite(surfWriter(), faces);
                 }
             }
         }
@@ -859,9 +874,9 @@ Foam::label Foam::checkGeometry
                 << " under-determined cells to set " << cells.name() << endl;
             cells.instance() = mesh.pointsInstance();
             cells.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), cells);
+                mergeAndWrite(surfWriter(), cells);
             }
         }
     }
@@ -879,9 +894,9 @@ Foam::label Foam::checkGeometry
                 << " concave cells to set " << cells.name() << endl;
             cells.instance() = mesh.pointsInstance();
             cells.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), cells);
+                mergeAndWrite(surfWriter(), cells);
             }
         }
     }
@@ -900,9 +915,9 @@ Foam::label Foam::checkGeometry
                 << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), faces);
+                mergeAndWrite(surfWriter(), faces);
             }
         }
     }
@@ -921,9 +936,9 @@ Foam::label Foam::checkGeometry
                 << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), faces);
+                mergeAndWrite(surfWriter(), faces);
             }
         }
     }
@@ -936,11 +951,16 @@ Foam::label Foam::checkGeometry
         const word procAndTime(Foam::name(Pstream::myProcNo()) + "_" + tmName);
 
         autoPtr<surfaceWriter> patchWriter;
-        if (!writer.valid())
+        if (!surfWriter.valid())
         {
             patchWriter.reset(new vtkSurfaceWriter());
         }
-        const surfaceWriter& wr = (writer.valid() ? writer() : patchWriter());
+        const surfaceWriter& wr =
+        (
+            surfWriter.valid()
+          ? surfWriter()
+          : patchWriter()
+        );
 
         forAll(pbm, patchi)
         {

applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H

@@ -1,6 +1,7 @@
 #include "label.H"
 #include "HashSet.H"
 #include "labelVector.H"
+#include "writer.H"
 
 namespace Foam
 {
@@ -26,6 +27,7 @@ namespace Foam
     (
         const polyMesh& mesh,
         const bool allGeometry,
-        const autoPtr<surfaceWriter>&
+        const autoPtr<surfaceWriter>&,
+        const autoPtr<writer<scalar>>&
     );
 }

applications/utilities/mesh/manipulation/checkMesh/checkMesh.C

@@ -47,7 +47,8 @@ Usage
       - \par -writeSets \<surfaceFormat\>
         Reconstruct all cellSets and faceSets geometry and write to
         postProcessing directory according to surfaceFormat
-        (e.g. vtk or ensight)
+        (e.g. vtk or ensight). Additionally reconstructs all pointSets and
+        writes as vtk format.
 
 \*---------------------------------------------------------------------------*/
 
@@ -56,7 +57,8 @@ Usage
 #include "Time.H"
 #include "polyMesh.H"
 #include "globalMeshData.H"
-#include "vtkSurfaceWriter.H"
+#include "surfaceWriter.H"
+#include "vtkSetWriter.H"
 
 #include "checkTools.H"
 #include "checkTopology.H"
@@ -94,7 +96,7 @@ int main(int argc, char *argv[])
     argList::addOption
     (
         "writeSets",
-        "<surfaceFormat>"
+        "surfaceFormat",
         "reconstruct and write all faceSets and cellSets in selected format"
     );
 
@@ -156,10 +158,12 @@ int main(int argc, char *argv[])
     }
 
 
-    autoPtr<surfaceWriter> writer;
+    autoPtr<surfaceWriter> surfWriter;
+    autoPtr<writer<scalar>> setWriter;
     if (writeSets)
     {
-        writer = surfaceWriter::New(surfaceFormat);
+        surfWriter = surfaceWriter::New(surfaceFormat);
+        setWriter = writer<scalar>::New(vtkSetWriter<scalar>::typeName);
     }
 
 
@@ -192,15 +196,22 @@ int main(int argc, char *argv[])
                     mesh,
                     allTopology,
                     allGeometry,
-                    writer
+                    surfWriter,
+                    setWriter
                 );
             }
 
-            nFailedChecks += checkGeometry(mesh, allGeometry, writer);
+            nFailedChecks += checkGeometry
+            (
+                mesh,
+                allGeometry,
+                surfWriter,
+                setWriter
+            );
 
             if (meshQuality)
             {
-                nFailedChecks += checkMeshQuality(mesh, qualDict(), writer);
+                nFailedChecks += checkMeshQuality(mesh, qualDict(), surfWriter);
             }
 
 
@@ -221,11 +232,17 @@ int main(int argc, char *argv[])
         {
             Info<< "Time = " << runTime.timeName() << nl << endl;
 
-            label nFailedChecks = checkGeometry(mesh, allGeometry, writer);
+            label nFailedChecks = checkGeometry
+            (
+                mesh,
+                allGeometry,
+                surfWriter,
+                setWriter
+            );
 
             if (meshQuality)
             {
-                nFailedChecks += checkMeshQuality(mesh, qualDict(), writer);
+                nFailedChecks += checkMeshQuality(mesh, qualDict(), surfWriter);
             }
 
 

applications/utilities/mesh/manipulation/checkMesh/checkTools.C

@@ -33,6 +33,7 @@ License
 #include "tetWedgeMatcher.H"
 #include "tetMatcher.H"
 #include "IOmanip.H"
+#include "pointSet.H"
 #include "faceSet.H"
 #include "cellSet.H"
 #include "Time.H"
@@ -369,6 +370,118 @@ void Foam::mergeAndWrite
       / set.name()
     );
 
-
     mergeAndWrite(mesh, writer, set.name(), setPatch, outputDir);
 }
+
+
+void Foam::mergeAndWrite
+(
+    const writer<scalar>& writer,
+    const pointSet& set
+)
+{
+    const polyMesh& mesh = refCast<const polyMesh>(set.db());
+
+    pointField mergedPts;
+    labelList mergedIDs;
+
+    if (Pstream::parRun())
+    {
+        // Note: we explicitly do not merge the points
+        // (mesh.globalData().mergePoints etc) since this might
+        // hide any synchronisation problem
+
+        globalIndex globalNumbering(mesh.nPoints());
+
+        mergedPts.setSize(returnReduce(set.size(), sumOp<label>()));
+        mergedIDs.setSize(mergedPts.size());
+
+        labelList setPointIDs(set.sortedToc());
+
+        // Get renumbered local data
+        pointField myPoints(mesh.points(), setPointIDs);
+        labelList myIDs(setPointIDs.size());
+        forAll(setPointIDs, i)
+        {
+            myIDs[i] = globalNumbering.toGlobal(setPointIDs[i]);
+        }
+
+        if (Pstream::master())
+        {
+            // Insert master data first
+            label pOffset = 0;
+            SubList<point>(mergedPts, myPoints.size(), pOffset) = myPoints;
+            SubList<label>(mergedIDs, myIDs.size(), pOffset) = myIDs;
+            pOffset += myPoints.size();
+
+            // Receive slave ones
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+
+                pointField slavePts(fromSlave);
+                labelList slaveIDs(fromSlave);
+
+                SubList<point>(mergedPts, slavePts.size(), pOffset) = slavePts;
+                SubList<label>(mergedIDs, slaveIDs.size(), pOffset) = slaveIDs;
+                pOffset += slaveIDs.size();
+            }
+        }
+        else
+        {
+            // Construct processor stream with estimate of size. Could
+            // be improved.
+            OPstream toMaster
+            (
+                Pstream::scheduled,
+                Pstream::masterNo(),
+                myPoints.byteSize() + myIDs.byteSize()
+            );
+            toMaster << myPoints << myIDs;
+        }
+    }
+    else
+    {
+        mergedIDs = set.sortedToc();
+        mergedPts = pointField(mesh.points(), mergedIDs);
+    }
+
+
+    // Write with scalar pointID
+    if (Pstream::master())
+    {
+        scalarField scalarPointIDs(mergedIDs.size());
+        forAll(mergedIDs, i)
+        {
+            scalarPointIDs[i] = 1.0*mergedIDs[i];
+        }
+
+        coordSet points(set.name(), "distance", mergedPts, mag(mergedPts));
+
+        List<const scalarField*> flds(1, &scalarPointIDs);
+
+        wordList fldNames(1, "pointID");
+
+        // Output e.g. pointSet p0 to
+        // postProcessing/<time>/p0.vtk
+        const fileName outputDir
+        (
+            set.time().path()
+          / (Pstream::parRun() ? ".." : "")
+          / "postProcessing"
+          / mesh.pointsInstance()
+          // set.name()
+        );
+        mkDir(outputDir);
+
+        fileName outputFile(outputDir/writer.getFileName(points, wordList()));
+        //fileName outputFile(outputDir/set.name());
+
+        OFstream os(outputFile);
+
+        writer.write(points, fldNames, flds, os);
+    }
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/manipulation/checkMesh/checkTools.H

@@ -1,10 +1,12 @@
 #include "scalar.H"
 #include "indirectPrimitivePatch.H"
+#include "writer.H"
 
 namespace Foam
 {
     class polyMesh;
     class surfaceWriter;
+    class pointSet;
     class faceSet;
     class cellSet;
     class fileName;
@@ -23,11 +25,18 @@ namespace Foam
         const fileName& outputDir
     );
 
-    //- Write vtk representation of (assembled) faceSet to vtk file in
+    //- Write vtk representation of (assembled) faceSet to surface file in
     //  postProcessing/ directory
     void mergeAndWrite(const surfaceWriter&, const faceSet&);
 
-    //- Write vtk representation of (assembled) cellSet to vtk file in
+    //- Write vtk representation of (assembled) cellSet to surface file in
     //  postProcessing/ directory
     void mergeAndWrite(const surfaceWriter&, const cellSet&);
+
+    //- Write vtk representation of (assembled) pointSet to 'set' file in
+    //  postProcessing/ directory
+    void mergeAndWrite(const writer<scalar>&, const pointSet&);
 }
+
+
+// ************************************************************************* //

applications/utilities/mesh/manipulation/checkMesh/checkTopology.C

@@ -43,7 +43,8 @@ Foam::label Foam::checkTopology
     const polyMesh& mesh,
     const bool allTopology,
     const bool allGeometry,
-    const autoPtr<surfaceWriter>& writer
+    const autoPtr<surfaceWriter>& surfWriter,
+    const autoPtr<writer<scalar>>& setWriter
 )
 {
     label noFailedChecks = 0;
@@ -129,9 +130,9 @@ Foam::label Foam::checkTopology
                 << " illegal cells to set " << cells.name() << endl;
             cells.instance() = mesh.pointsInstance();
             cells.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), cells);
+                mergeAndWrite(surfWriter(), cells);
             }
 
         }
@@ -154,6 +155,10 @@ Foam::label Foam::checkTopology
                 << " unused points to set " << points.name() << endl;
             points.instance() = mesh.pointsInstance();
             points.write();
+            if (setWriter.valid())
+            {
+                mergeAndWrite(setWriter, points);
+            }
         }
     }
 
@@ -172,9 +177,9 @@ Foam::label Foam::checkTopology
                 << " unordered faces to set " << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), faces);
+                mergeAndWrite(surfWriter(), faces);
             }
         }
     }
@@ -192,9 +197,9 @@ Foam::label Foam::checkTopology
                 << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), faces);
+                mergeAndWrite(surfWriter(), faces);
             }
         }
     }
@@ -213,9 +218,9 @@ Foam::label Foam::checkTopology
                 << endl;
             cells.instance() = mesh.pointsInstance();
             cells.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), cells);
+                mergeAndWrite(surfWriter(), cells);
             }
 
         }
@@ -237,9 +242,9 @@ Foam::label Foam::checkTopology
                 << faces.name() << endl;
             faces.instance() = mesh.pointsInstance();
             faces.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), faces);
+                mergeAndWrite(surfWriter(), faces);
             }
         }
     }
@@ -292,9 +297,9 @@ Foam::label Foam::checkTopology
                 << endl;
             oneCells.instance() = mesh.pointsInstance();
             oneCells.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), oneCells);
+                mergeAndWrite(surfWriter(), oneCells);
             }
         }
 
@@ -308,9 +313,9 @@ Foam::label Foam::checkTopology
                 << endl;
             twoCells.instance() = mesh.pointsInstance();
             twoCells.write();
-            if (writer.valid())
+            if (surfWriter.valid())
             {
-                mergeAndWrite(writer(), twoCells);
+                mergeAndWrite(surfWriter(), twoCells);
             }
         }
     }
@@ -349,6 +354,13 @@ Foam::label Foam::checkTopology
             ctr.write();
 
 
+            // Points in multiple regions
+            pointSet points
+            (
+                mesh,
+                "multiRegionPoints",
+                mesh.nPoints()/1000
+            );
 
             // Is region disconnected
             boolList regionDisconnected(rs.nRegions(), true);
@@ -386,6 +398,7 @@ Foam::label Foam::checkTopology
                             regionDisconnected[regionI] = false;
                             regionDisconnected[pRegion] = false;
                             pRegion = -2;
+                            points.insert(f[fp]);
                         }
                     }
                 }
@@ -437,6 +450,19 @@ Foam::label Foam::checkTopology
 
                 cellRegions[i].write();
             }
+
+            label nPoints = returnReduce(points.size(), sumOp<label>());
+            if (nPoints)
+            {
+                Info<< "  <<Writing " << nPoints
+                    << " points that are in multiple regions to set "
+                    << points.name() << endl;
+                points.write();
+                if (setWriter.valid())
+                {
+                    mergeAndWrite(setWriter, points);
+                }
+            }
         }
     }
 

applications/utilities/mesh/manipulation/checkMesh/checkTopology.H

@@ -1,5 +1,6 @@
 #include "label.H"
 #include "autoPtr.H"
+#include "writer.H"
 
 namespace Foam
 {
@@ -11,6 +12,7 @@ namespace Foam
         const polyMesh&,
         const bool,
         const bool,
-        const autoPtr<surfaceWriter>&
+        const autoPtr<surfaceWriter>&,
+        const autoPtr<writer<scalar>>&
     );
 }