checkMesh, functionObjects::checkMesh: Updated and made consistent

Now both the checkMesh utility and functionObject operate in the same manner
with the same controls, executing the same checkGeometry and checkTopology
functions from the new meshCheck library.  The controls have been updated and
made more consistent and flexible, in particular by the addition of optional
user specification for the non-orthogonality and skewness error thresholds:

Application
    checkMesh

Description
    Checks validity of a mesh.

Usage
    \b checkMesh [OPTION]

    Options:
      - \par noTopology
        Skip checking the mesh topology

      - \par -allTopology
        Check all (including non finite-volume specific) addressing

      - \par -allGeometry
        Check all (including non finite-volume specific) geometry

      - \par -meshQuality
        Check against user defined (in \a system/meshQualityDict) quality
        settings

      - \par -region \<name\>
        Specify an alternative mesh region.

      - \par -writeSurfaces
        Reconstruct cellSets and faceSets of problem faces and write to
        postProcessing directory.

      - \par -surfaceFormat <format>
        Format used to write the cellSets and faceSets surfaces
        e.g. vtk or ensight.

      - \par -writeSets
        Reconstruct pointSets of problem points nd write to
        postProcessing directory.

      - \par -setFormat <format>
        Format used to write the pointSets
        e.g. vtk or ensight.

      - \par -nonOrthThreshold <threshold value in degrees>
        Threshold in degrees for reporting non-orthogonality errors,
        default: 70"

      - \par -skewThreshold <threshold value>
        Threshold for reporting skewness errors, default: 4.

Class
    Foam::functionObjects::checkMesh

Description
    Executes primitiveMesh::checkMesh(true) every execute time for which the
    mesh changed, i.e. moved or changed topology.

    Useful to check the correctness of changing and morphing meshes.

Usage
    \table
        Property      | Description               | Required     | Default value
        type          | type name: checkMesh      | yes          |
        noTopology    | Skip checking the mesh topology | no     | false
        allTopology   | Check all addressing      | no           | false
        allGeometry   | Check all geometry        | no           | false
        writeSurfaces | Reconstruct and write problem faces | no | false
        surfaceFormat | Format for problem faceSets | no         | vtk
        writeSets     | Reconstruct and write problem points | no | false
        setFormat     | Format used to write the problem pointSets | no | vtk
        nonOrthThreshold | Threshold for non-orthogonality errors | no | 70 deg
        skewThreshold | Threshold for reporting skewness errors | no | 4
    \endtable

    Example of checkMesh specification:
    \verbatim
    checkMesh
    {
        type            checkMesh;
        libs            ("libutilityFunctionObjects.so");

        executeControl  timeStep;
        executeInterval 10;

        allGeometry     true;
        allTopology     true;

        writeSurfaces   true;
        surfaceFormat   vtk;

        writeSets       true;
        setFormat       vtk;
    }
    \endverbatim

    or using the standard configuration file:

    \verbatim
    #includeFunc checkMesh(executeInterval=10, allGeometry=true)
    \endverbatim

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

@@ -31,24 +31,44 @@ Usage
     \b checkMesh [OPTION]
 
     Options:
-      - \par -allGeometry
+      - \par noTopology
-        Checks all (including non finite-volume specific) geometry
+        Skip checking the mesh topology
 
       - \par -allTopology
-        Checks all (including non finite-volume specific) addressing
+        Check all (including non finite-volume specific) addressing
+
+      - \par -allGeometry
+        Check all (including non finite-volume specific) geometry
 
       - \par -meshQuality
-        Checks against user defined (in \a system/meshQualityDict) quality
+        Check against user defined (in \a system/meshQualityDict) quality
         settings
 
       - \par -region \<name\>
         Specify an alternative mesh region.
 
-      - \par -writeSets \<surfaceFormat\>
+      - \par -writeSurfaces
-        Reconstruct all cellSets and faceSets geometry and write to
+        Reconstruct cellSets and faceSets of problem faces and write to
-        postProcessing directory according to surfaceFormat
+        postProcessing directory.
-        (e.g. vtk or ensight). Additionally reconstructs all pointSets and
+
-        writes as vtk format.
+      - \par -surfaceFormat <format>
+        Format used to write the cellSets and faceSets surfaces
+        e.g. vtk or ensight.
+
+      - \par -writeSets
+        Reconstruct pointSets of problem points nd write to
+        postProcessing directory.
+
+      - \par -setFormat <format>
+        Format used to write the pointSets
+        e.g. vtk or ensight.
+
+      - \par -nonOrthThreshold <threshold value in degrees>
+        Threshold in degrees for reporting non-orthogonality errors,
+        default: 70"
+
+      - \par -skewThreshold <threshold value>
+        Threshold for reporting skewness errors, default: 4.
 
 \*---------------------------------------------------------------------------*/
 
@@ -57,7 +77,7 @@ Usage
 #include "Time.H"
 #include "polyMesh.H"
 #include "globalMeshData.H"
-#include "surfaceWriter.H"
+#include "vtkSurfaceWriter.H"
 #include "vtkSetWriter.H"
 #include "IOdictionary.H"
 
@@ -92,11 +112,40 @@ int main(int argc, char *argv[])
         "meshQuality",
         "read user-defined mesh quality criterions from system/meshQualityDict"
     );
+    argList::addBoolOption
+    (
+        "writeSurfaces",
+        "reconstruct and write faceSets and cellSets of the problem faces"
+    );
     argList::addOption
     (
-        "writeSets",
         "surfaceFormat",
-        "reconstruct and write all faceSets and cellSets in selected format"
+        "surfaceFormat",
+        "Format for faceSets and cellSets of the problem faces, defaults to vtk"
+    );
+    argList::addBoolOption
+    (
+        "writeSets",
+        "reconstruct and write pointSets of the problem points"
+    );
+    argList::addOption
+    (
+        "setFormat",
+        "setFormat",
+        "Format for pointSets of the problem points, defaults to vtk"
+    );
+    argList::addOption
+    (
+        "nonOrthThreshold",
+        "nonOrthThreshold",
+        "Threshold in degrees "
+        "for reporting non-orthogonality errors, default: 70"
+    );
+    argList::addOption
+    (
+        "skewThreshold",
+        "skewThreshold",
+        "Threshold for reporting non-orthogonality errors, default: 4"
     );
 
     #include "setRootCase.H"
@@ -108,9 +157,20 @@ int main(int argc, char *argv[])
     const bool allGeometry = args.optionFound("allGeometry");
     const bool allTopology = args.optionFound("allTopology");
     const bool meshQuality = args.optionFound("meshQuality");
+    const bool writeSurfaces = args.optionFound("writeSurfaces");
+    const bool writeSets = args.optionFound("writeSets");
 
-    word surfaceFormat;
+    word surfaceFormat(vtkSurfaceWriter::typeName);
-    const bool writeSets = args.optionReadIfPresent("writeSets", surfaceFormat);
+    args.optionReadIfPresent("surfaceFormat", surfaceFormat);
+
+    word setFormat(vtkSetWriter::typeName);
+    args.optionReadIfPresent("surfaceFormat", setFormat);
+
+    scalar nonOrthThreshold = 70;
+    args.optionReadIfPresent("nonOrthThreshold", nonOrthThreshold);
+
+    scalar skewThreshold = 4;
+    args.optionReadIfPresent("skewThreshold", skewThreshold);
 
     if (noTopology)
     {
@@ -129,11 +189,17 @@ int main(int argc, char *argv[])
     {
         Info<< "Enabling user-defined geometry checks." << nl << endl;
     }
+    if (writeSurfaces)
+    {
+        Info<< "Reconstructing and writing surface representation of the "
+            << "faceSets and cellSets of problem faces in "
+            << surfaceFormat << " format" << nl << endl;
+    }
     if (writeSets)
     {
-        Info<< "Reconstructing and writing " << surfaceFormat
+        Info<< "Reconstructing and writing the problem points in "
-            << " representation"
+            << setFormat << " format"
-            << " of all faceSets and cellSets." << nl << endl;
+            << nl << endl;
     }
 
 
@@ -156,20 +222,23 @@ int main(int argc, char *argv[])
         );
     }
 
-
+    autoPtr<Foam::surfaceWriter> surfaceWriter;
-    autoPtr<surfaceWriter> surfWriter;
+    if (writeSurfaces)
-    autoPtr<Foam::setWriter> setWriter;
-    if (writeSets)
     {
-        surfWriter = surfaceWriter::New
+        surfaceWriter = surfaceWriter::New
         (
             surfaceFormat,
             mesh.time().writeFormat(),
             mesh.time().writeCompression()
         );
+    }
+
+    autoPtr<Foam::setWriter> setWriter;
+    if (writeSets)
+    {
         setWriter = Foam::setWriter::New
         (
-            vtkSetWriter::typeName,
+            setFormat,
             mesh.time().writeFormat(),
             mesh.time().writeCompression()
         );
@@ -204,7 +273,7 @@ int main(int argc, char *argv[])
                 (
                     mesh,
                     allTopology,
-                    surfWriter,
+                    surfaceWriter,
                     setWriter
                 );
             }
@@ -213,13 +282,16 @@ int main(int argc, char *argv[])
             (
                 mesh,
                 allGeometry,
-                surfWriter,
+                nonOrthThreshold,
+                skewThreshold,
+                surfaceWriter,
                 setWriter
             );
 
             if (meshQuality)
             {
-                nFailedChecks += checkMeshQuality(mesh, qualDict(), surfWriter);
+                nFailedChecks +=
+                    checkMeshQuality(mesh, qualDict(), surfaceWriter);
             }
 
 
@@ -244,13 +316,16 @@ int main(int argc, char *argv[])
             (
                 mesh,
                 allGeometry,
-                surfWriter,
+                nonOrthThreshold,
+                skewThreshold,
+                surfaceWriter,
                 setWriter
             );
 
             if (meshQuality)
             {
-                nFailedChecks += checkMeshQuality(mesh, qualDict(), surfWriter);
+                nFailedChecks +=
+                    checkMeshQuality(mesh, qualDict(), surfaceWriter);
             }
 
 

src/dynamicMesh/polyTopoChange/attachPolyTopoChanger/attachPolyTopoChanger.C

@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "attachPolyTopoChanger.H"
-#include "meshCheck.H"
+#include "polyMesh.H"
 #include "polyTopoChange.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -119,8 +119,6 @@ void Foam::attachPolyTopoChanger::attach(const bool removeEmptyPatches)
         Pout<< "void attachPolyTopoChanger::attach(): "
             << "Finished attaching mesh" << endl;
     }
-
-    meshCheck::checkMesh(mesh_);
 }
 
 

src/functionObjects/utilities/Make/options

@@ -1,7 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshCheck/lnInclude
+    -I$(LIB_SRC)/meshCheck/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
-    -lmeshCheck
+    -lmeshCheck \
+    -lsampling

src/functionObjects/utilities/checkMesh/checkMesh.C

@@ -26,6 +26,8 @@ License
 #include "checkMesh.H"
 #include "fvMesh.H"
 #include "meshCheck.H"
+#include "vtkSetWriter.H"
+#include "vtkSurfaceWriter.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -63,11 +65,99 @@ Foam::functionObjects::checkMesh::~checkMesh()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::functionObjects::checkMesh::read(const dictionary& dict)
+{
+    noTopology_ =  dict.lookupOrDefault("noTopology", false);
+    allGeometry_ = dict.lookupOrDefault("allGeometry", false);
+    allTopology_ = dict.lookupOrDefault("allTopology", false);
+
+    writeSurfaces_ = dict.lookupOrDefault("writeSurfaces", false);
+    if (writeSurfaces_)
+    {
+        surfaceFormat_ = dict.lookupOrDefault<word>
+        (
+            "surfaceFormat",
+            vtkSurfaceWriter::typeName
+        );
+    }
+
+    writeSets_ = dict.lookupOrDefault("writeSets", false);
+    if (writeSets_)
+    {
+        setFormat_ = dict.lookupOrDefault<word>
+        (
+            "setFormat",
+            vtkSetWriter::typeName
+        );
+    }
+
+    nonOrthThreshold_ = dict.lookupOrDefault("nonOrthThreshold", 70.0);
+    skewThreshold_ = dict.lookupOrDefault("skewThreshold", 4.0);
+
+    return functionObject::read(dict);
+}
+
+
 bool Foam::functionObjects::checkMesh::execute()
 {
     if (mesh_.changing())
     {
-        return meshCheck::checkMesh(mesh_, true);
+        autoPtr<surfaceWriter> surfWriter;
+        if (writeSurfaces_)
+        {
+            surfWriter = surfaceWriter::New
+            (
+                surfaceFormat_,
+                mesh_.time().writeFormat(),
+                mesh_.time().writeCompression()
+            );
+        }
+
+        autoPtr<setWriter> setWriter;
+        if (writeSets_)
+        {
+            setWriter = setWriter::New
+            (
+                setFormat_,
+                mesh_.time().writeFormat(),
+                mesh_.time().writeCompression()
+            );
+        }
+
+        label nFailedChecks = 0;
+
+        if (!noTopology_)
+        {
+            nFailedChecks += meshCheck::checkTopology
+            (
+                mesh_,
+                allTopology_,
+                surfWriter,
+                setWriter
+            );
+        }
+
+        nFailedChecks += meshCheck::checkGeometry
+        (
+            mesh_,
+            allGeometry_,
+            nonOrthThreshold_,
+            skewThreshold_,
+            surfWriter,
+            setWriter
+        );
+
+        if (nFailedChecks == 0)
+        {
+            Info<< "\n    Mesh OK.\n" << endl;
+        }
+        else
+        {
+            Info<< "\n    Failed " << nFailedChecks << " mesh checks.\n"
+                << endl;
+        }
+
+        return nFailedChecks == 0;
     }
     else
     {

src/functionObjects/utilities/checkMesh/checkMesh.H

@@ -30,6 +30,21 @@ Description
 
     Useful to check the correctness of changing and morphing meshes.
 
+Usage
+    \table
+        Property      | Description               | Required     | Default value
+        type          | type name: checkMesh      | yes          |
+        noTopology    | Skip checking the mesh topology | no     | false
+        allTopology   | Check all addressing      | no           | false
+        allGeometry   | Check all geometry        | no           | false
+        writeSurfaces | Reconstruct and write problem faces | no | false
+        surfaceFormat | Format for problem faceSets | no         | vtk
+        writeSets     | Reconstruct and write problem points | no | false
+        setFormat     | Format used to write the problem pointSets | no | vtk
+        nonOrthThreshold | Threshold for non-orthogonality errors | no | 70 deg
+        skewThreshold | Threshold for reporting skewness errors | no | 4
+    \endtable
+
     Example of checkMesh specification:
     \verbatim
     checkMesh
@@ -39,11 +54,22 @@ Description
 
         executeControl  timeStep;
         executeInterval 10;
+
+        allGeometry     true;
+        allTopology     true;
+
+        writeSurfaces   true;
+        surfaceFormat   vtk;
+
+        writeSets       true;
+        setFormat       vtk;
     }
     \endverbatim
+
     or using the standard configuration file:
+
     \verbatim
-    #includeFunc checkMesh(executeInterval=10)
+    #includeFunc checkMesh(executeInterval=10, allGeometry=true)
     \endverbatim
 
 SourceFiles
@@ -73,6 +99,20 @@ class checkMesh
 {
     // Private Data
 
+        bool noTopology_;
+        bool allGeometry_;
+        bool allTopology_;
+
+        bool writeSurfaces_;
+        word surfaceFormat_;
+
+        bool writeSets_;
+        word setFormat_;
+
+        scalar nonOrthThreshold_;
+        scalar skewThreshold_;
+
+
 public:
 
     //- Runtime type information
@@ -96,6 +136,9 @@ public:
 
     // Member Functions
 
+        //- Read the checkMesh controls
+        virtual bool read(const dictionary&);
+
         //- Return the list of fields required
         virtual wordList fields() const
         {

src/meshCheck/checkGeometry.C

@@ -509,6 +509,8 @@ Foam::label Foam::meshCheck::checkGeometry
 (
     const polyMesh& mesh,
     const bool allGeometry,
+    const scalar nonOrthThreshold,
+    const scalar skewThreshold,
     const autoPtr<surfaceWriter>& surfWriter,
     const autoPtr<Foam::setWriter>& setWriter
 )
@@ -710,7 +712,16 @@ Foam::label Foam::meshCheck::checkGeometry
 
     {
         faceSet faces(mesh, "nonOrthoFaces", mesh.nFaces()/100+1);
-        if (meshCheck::checkFaceOrthogonality(mesh, true, &faces))
+        if
+        (
+            meshCheck::checkFaceOrthogonality
+            (
+                mesh,
+                nonOrthThreshold,
+                true,
+                &faces
+            )
+        )
         {
             noFailedChecks++;
         }
@@ -755,7 +766,16 @@ Foam::label Foam::meshCheck::checkGeometry
 
     {
         faceSet faces(mesh, "skewFaces", mesh.nFaces()/100+1);
-        if (meshCheck::checkFaceSkewness(mesh, true, &faces))
+        if
+        (
+            meshCheck::checkFaceSkewness
+            (
+                mesh,
+                skewThreshold,
+                true,
+                &faces
+            )
+        )
         {
             noFailedChecks++;
 

src/meshCheck/meshCheck.H

@@ -81,6 +81,8 @@ namespace meshCheck
     (
         const polyMesh& mesh,
         const bool allGeometry,
+        const scalar nonOrthThreshold,
+        const scalar skewThreshold,
         const autoPtr<surfaceWriter>&,
         const autoPtr<setWriter>&
     );

src/meshCheck/polyMeshCheck/polyMeshCheck.C

@@ -281,6 +281,7 @@ Foam::tmp<Foam::scalarField> Foam::meshCheck::volRatio
 bool Foam::meshCheck::checkFaceOrthogonality
 (
     const polyMesh& mesh,
+    const scalar nonOrthThreshold,
     const bool report,
     labelHashSet* setPtr
 )
@@ -305,7 +306,7 @@ bool Foam::meshCheck::checkFaceOrthogonality
 
     // Severe nonorthogonality threshold
     const scalar severeNonorthogonalityThreshold =
-        ::cos(degToRad(meshCheck::nonOrthThreshold));
+        ::cos(degToRad(nonOrthThreshold));
 
 
     scalar minDDotS = great;
@@ -374,7 +375,7 @@ bool Foam::meshCheck::checkFaceOrthogonality
         if (severeNonOrth > 0)
         {
             Info<< "   *Number of severely non-orthogonal (> "
-                << meshCheck::nonOrthThreshold << " degrees) faces: "
+                << nonOrthThreshold << " degrees) faces: "
                 << severeNonOrth << "." << endl;
         }
     }
@@ -404,6 +405,7 @@ bool Foam::meshCheck::checkFaceOrthogonality
 bool Foam::meshCheck::checkFaceSkewness
 (
     const polyMesh& mesh,
+    const scalar skewThreshold,
     const bool report,
     labelHashSet* setPtr
 )
@@ -441,7 +443,7 @@ bool Foam::meshCheck::checkFaceSkewness
     {
         // Check if the skewness vector is greater than the PN vector.
         // This does not cause trouble but is a good indication of a poor mesh.
-        if (skew[facei] > meshCheck::skewThreshold)
+        if (skew[facei] > skewThreshold)
         {
             if (setPtr)
             {

src/meshCheck/polyMeshCheck/polyMeshCheck.H

@@ -87,6 +87,7 @@ namespace meshCheck
     bool checkFaceOrthogonality
     (
         const polyMesh& mesh,
+        const scalar nonOrthThreshold,
         const bool report = false,
         labelHashSet* setPtr = nullptr
     );
@@ -95,6 +96,7 @@ namespace meshCheck
     bool checkFaceSkewness
     (
         const polyMesh& mesh,
+        const scalar skewThreshold,
         const bool report = false,
         labelHashSet* setPtr = nullptr
     );

src/meshCheck/polyMeshCheck/polyMeshCheckQuality.C

@@ -42,7 +42,7 @@ namespace meshCheck
 
 scalar checkNonOrtho
 (
-    const polyMesh& mesh,
+    const primitiveMesh& mesh,
     const bool report,
     const scalar severeNonorthogonalityThreshold,
     const label facei,
@@ -109,7 +109,7 @@ scalar checkNonOrtho
 
 bool checkFaceTet
 (
-    const polyMesh& mesh,
+    const primitiveMesh& mesh,
     const bool report,
     const scalar minTetQuality,
     const pointField& p,
@@ -161,7 +161,7 @@ bool checkFaceTet
 
 labelList getAffectedCells
 (
-    const polyMesh& mesh,
+    const primitiveMesh& mesh,
     const labelList& changedFaces
 )
 {
@@ -184,6 +184,7 @@ labelList getAffectedCells
     return affectedCells.toc();
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace meshCheck
@@ -1443,53 +1444,8 @@ bool Foam::meshCheck::checkFaceTwist
 
     const faceList& fcs = mesh.faces();
 
-
     label nWarped = 0;
 
-//    forAll(checkFaces, i)
-//    {
-//        label facei = checkFaces[i];
-//
-//        const face& f = fcs[facei];
-//
-//        scalar magArea = mag(faceAreas[facei]);
-//
-//        if (f.size() > 3 && magArea > vSmall)
-//        {
-//            const vector nf = faceAreas[facei] / magArea;
-//
-//            const point& fc = faceCentres[facei];
-//
-//            forAll(f, fpI)
-//            {
-//                vector triArea
-//                (
-//                    triPointRef
-//                    (
-//                        p[f[fpI]],
-//                        p[f.nextLabel(fpI)],
-//                        fc
-//                    ).area()
-//                );
-//
-//                scalar magTri = mag(triArea);
-//
-//                if (magTri > vSmall && ((nf & triArea/magTri) < minTwist))
-//                {
-//                    nWarped++;
-//
-//                    if (setPtr)
-//                    {
-//                        setPtr->insert(facei);
-//                    }
-//
-//                    break;
-//                }
-//            }
-//        }
-//    }
-
-
     const labelList& own = mesh.faceOwner();
     const labelList& nei = mesh.faceNeighbour();
     const polyBoundaryMesh& patches = mesh.boundaryMesh();

src/meshCheck/primitiveMeshCheck/primitiveMeshCheck.C

@@ -35,8 +35,6 @@ License
 
 Foam::scalar Foam::meshCheck::closedThreshold  = 1.0e-6;
 Foam::scalar Foam::meshCheck::aspectThreshold  = 1000;
-Foam::scalar Foam::meshCheck::nonOrthThreshold = 70;    // deg
-Foam::scalar Foam::meshCheck::skewThreshold    = 4;
 Foam::scalar Foam::meshCheck::planarCosAngle   = 1.0e-6;