From f276366a050e99812c5f4fa7e3d0ef13be90c511 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 21 Sep 2022 11:18:57 +0100
Subject: [PATCH] ENH: snappyHexMesh: add mesh-quality control for edge lengths
---
etc/caseDicts/meshQualityDict | 4 +-
.../motionSmoother/motionSmootherAlgoCheck.C | 105 +++++++++++++++++-
.../polyMeshGeometry/polyMeshGeometry.C | 88 +++++++++++++++
.../polyMeshGeometry/polyMeshGeometry.H | 20 +++-
4 files changed, 212 insertions(+), 5 deletions(-)
diff --git a/etc/caseDicts/meshQualityDict b/etc/caseDicts/meshQualityDict
index ada9935597e..b1e007fe1c0 100644
--- a/etc/caseDicts/meshQualityDict
+++ b/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;
+
// ************************************************************************* //
diff --git a/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C b/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
index 85c7a72d81f..1b3a8701021 100644
--- a/src/dynamicMesh/motionSmoother/motionSmootherAlgoCheck.C
+++ b/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;
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
index 77e655392ab..8113c2262dc 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C
+++ b/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
+ );
+}
+
+
// ************************************************************************* //
diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H
index 242c925e090..e8115e6c9d5 100644
--- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.H
+++ b/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;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
--
GitLab