ENH: snappyHexMesh: add mesh-quality control for edge lengths

etc/caseDicts/meshQualityDict

@@ -67,12 +67,14 @@ minVolRatio 0.01;
 //  for Fluent compatibility
 minTriangleTwist -1;
 
-
 //- If >0 : preserve cells with all points on the surface if the
 //  resulting volume after snapping (by approximation) is larger than
 //  minVolCollapseRatio times old volume (i.e. not collapsed to flat cell).
 //  If <0 : delete always.
 //minVolCollapseRatio 0.1;
 
+//- Minimum edge length. Set to <0 to disable
+minEdgeLength 0.001;
+
 
 // ************************************************************************* //

src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2014 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020,2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ License
 #include "motionSmootherAlgo.H"
 #include "polyMeshGeometry.H"
 #include "IOmanip.H"
+#include "pointSet.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
@@ -127,6 +128,13 @@ bool Foam::motionSmootherAlgo::checkMesh
     (
         get<scalar>(dict, "minDeterminant", dryRun, keyType::REGEX_RECURSIVE)
     );
+    const scalar minEdgeLength
+    (
+        dict.getOrDefault<scalar>
+        (
+            "minEdgeLength", -1, keyType::REGEX_RECURSIVE
+        )
+    );
 
 
     if (dryRun)
@@ -452,6 +460,40 @@ bool Foam::motionSmootherAlgo::checkMesh
         nWrongFaces = nNewWrongFaces;
     }
 
+    if (minEdgeLength >= 0)
+    {
+        pointSet edgePoints(mesh, "smallEdgePoints", mesh.nPoints()/1000);
+        polyMeshGeometry::checkEdgeLength
+        (
+            report,
+            minEdgeLength,
+            mesh,
+            checkFaces,
+            &edgePoints
+        );
+
+        const auto& pointFaces = mesh.pointFaces();
+
+        label nNewWrongFaces = 0;
+        for (const label pointi : edgePoints)
+        {
+            const auto& pFaces = pointFaces[pointi];
+            for (const label facei : pFaces)
+            {
+                if (wrongFaces.insert(facei))
+                {
+                    nNewWrongFaces++;
+                }
+            }
+        }
+
+        Info<< "    faces with edge length < "
+            << setw(5) << minEdgeLength << "                         : "
+            << returnReduce(nNewWrongFaces, sumOp<label>()) << endl;
+
+        nWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
+    }
+
     //Pout.setf(ios_base::right);
 
     return nWrongFaces > 0;
@@ -539,11 +581,23 @@ bool Foam::motionSmootherAlgo::checkMesh
     );
     const scalar maxIntSkew
     (
-        get<scalar>(dict, "maxInternalSkewness", dryRun, keyType::REGEX_RECURSIVE)
+        get<scalar>
+        (
+            dict,
+            "maxInternalSkewness",
+            dryRun,
+            keyType::REGEX_RECURSIVE
+        )
     );
     const scalar maxBounSkew
     (
-        get<scalar>(dict, "maxBoundarySkewness", dryRun, keyType::REGEX_RECURSIVE)
+        get<scalar>
+        (
+            dict,
+            "maxBoundarySkewness",
+            dryRun,
+            keyType::REGEX_RECURSIVE
+        )
     );
     const scalar minWeight
     (
@@ -572,6 +626,13 @@ bool Foam::motionSmootherAlgo::checkMesh
     (
         get<scalar>(dict, "minDeterminant", dryRun, keyType::REGEX_RECURSIVE)
     );
+    const scalar minEdgeLength
+    (
+        dict.getOrDefault<scalar>
+        (
+            "minEdgeLength", -1, keyType::REGEX_RECURSIVE
+        )
+    );
 
     if (dryRun)
     {
@@ -865,6 +926,44 @@ bool Foam::motionSmootherAlgo::checkMesh
         nWrongFaces = nNewWrongFaces;
     }
 
+    if (minEdgeLength >= 0)
+    {
+        pointSet edgePoints
+        (
+            meshGeom.mesh(),
+            "smallEdgePoints",
+            meshGeom.mesh().nPoints()/1000
+        );
+        meshGeom.checkEdgeLength
+        (
+            report,
+            minEdgeLength,
+            checkFaces,
+            &edgePoints
+        );
+
+        const auto& pointFaces = meshGeom.mesh().pointFaces();
+
+        label nNewWrongFaces = 0;
+        for (const label pointi : edgePoints)
+        {
+            const auto& pFaces = pointFaces[pointi];
+            for (const label facei : pFaces)
+            {
+                if (wrongFaces.insert(facei))
+                {
+                    nNewWrongFaces++;
+                }
+            }
+        }
+
+        Info<< "    faces with edge length < "
+            << setw(5) << minEdgeLength << "                         : "
+            << returnReduce(nNewWrongFaces, sumOp<label>()) << endl;
+
+        nWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());
+    }
+
     //Pout.setf(ios_base::right);
 
     return nWrongFaces > 0;

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C

@@ -2085,6 +2085,75 @@ bool Foam::polyMeshGeometry::checkCellDeterminant
 }
 
 
+bool Foam::polyMeshGeometry::checkEdgeLength
+(
+    const bool report,
+    const scalar minEdgeLength,
+    const polyMesh& mesh,
+    const labelList& checkFaces,
+    labelHashSet* setPtr
+)
+{
+    const scalar reportLenSqr(Foam::sqr(minEdgeLength));
+
+    const pointField& points = mesh.points();
+    const faceList& faces = mesh.faces();
+
+    scalar minLenSqr = sqr(GREAT);
+    scalar maxLenSqr = -sqr(GREAT);
+
+    label nSmall = 0;
+
+    for (const label facei : checkFaces)
+    {
+        const face& f = faces[facei];
+
+        forAll(f, fp)
+        {
+            label fp1 = f.fcIndex(fp);
+
+            scalar magSqrE = magSqr(points[f[fp]] - points[f[fp1]]);
+
+            if (setPtr && magSqrE < reportLenSqr)
+            {
+                if (setPtr->insert(f[fp]) || setPtr->insert(f[fp1]))
+                {
+                    nSmall++;
+                }
+            }
+
+            minLenSqr = min(minLenSqr, magSqrE);
+            maxLenSqr = max(maxLenSqr, magSqrE);
+        }
+    }
+
+    reduce(minLenSqr, minOp<scalar>());
+    reduce(maxLenSqr, maxOp<scalar>());
+
+    reduce(nSmall, sumOp<label>());
+
+    if (nSmall > 0)
+    {
+        if (report)
+        {
+            Info<< "   *Edges too small, min/max edge length = "
+                << sqrt(minLenSqr) << " " << sqrt(maxLenSqr)
+                << ", number too small: " << nSmall << endl;
+        }
+
+        return true;
+    }
+
+    if (report)
+    {
+        Info<< "    Min/max edge length = "
+            << sqrt(minLenSqr) << " " << sqrt(maxLenSqr) << " OK." << endl;
+    }
+
+    return false;
+}
+
+
 bool Foam::polyMeshGeometry::checkFaceDotProduct
 (
     const bool report,
@@ -2363,4 +2432,23 @@ bool Foam::polyMeshGeometry::checkCellDeterminant
 }
 
 
+bool Foam::polyMeshGeometry::checkEdgeLength
+(
+    const bool report,
+    const scalar minEdgeLength,
+    const labelList& checkFaces,
+    labelHashSet* setPtr
+) const
+{
+    return checkEdgeLength
+    (
+        report,
+        minEdgeLength,
+        mesh_,
+        checkFaces,
+        setPtr
+    );
+}
+
+
 // ************************************************************************* //

src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020 OpenCFD Ltd.
+    Copyright (C) 2020,2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -343,6 +343,16 @@ public:
                 labelHashSet* setPtr
             );
 
+            //- Small edges. Optionally collects points of small edges
+            static bool checkEdgeLength
+            (
+                const bool report,
+                const scalar minEdgeLength,
+                const polyMesh& mesh,
+                const labelList& checkFaces,
+                labelHashSet* pointSetPtr
+            );
+
 
         // Checking of selected faces with local geometry. Uses above static
         // functions. Parallel aware.
@@ -446,6 +456,14 @@ public:
                 const labelList& affectedCells,
                 labelHashSet* setPtr
             ) const;
+
+            bool checkEdgeLength
+            (
+                const bool report,
+                const scalar minEdgeLength,
+                const labelList& checkFaces,
+                labelHashSet* pointSetPtr
+            ) const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //