From 4ab99793058a1facf3ea49f99097e868e14c1517 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Tue, 24 Oct 2023 14:36:35 +0100 Subject: [PATCH] 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 \ Specify an alternative mesh region. - \par -writeSurfaces Reconstruct cellSets and faceSets of problem faces and write to postProcessing directory. - \par -surfaceFormat 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 used to write the pointSets e.g. vtk or ensight. - \par -nonOrthThreshold Threshold in degrees for reporting non-orthogonality errors, default: 70" - \par -skewThreshold 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 --- .../mesh/manipulation/checkMesh/checkMesh.C | 131 ++++++++++++++---- .../attachPolyTopoChanger.C | 4 +- src/functionObjects/utilities/Make/options | 6 +- .../utilities/checkMesh/checkMesh.C | 92 +++++++++++- .../utilities/checkMesh/checkMesh.H | 45 +++++- src/meshCheck/checkGeometry.C | 24 +++- src/meshCheck/meshCheck.H | 2 + src/meshCheck/polyMeshCheck/polyMeshCheck.C | 8 +- src/meshCheck/polyMeshCheck/polyMeshCheck.H | 2 + .../polyMeshCheck/polyMeshCheckQuality.C | 52 +------ .../primitiveMeshCheck/primitiveMeshCheck.C | 2 - 11 files changed, 278 insertions(+), 90 deletions(-) diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C index 64de1960a8..9e9f24f5b6 100644 --- a/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C +++ b/applications/utilities/mesh/manipulation/checkMesh/checkMesh.C @@ -31,24 +31,44 @@ Usage \b checkMesh [OPTION] Options: - - \par -allGeometry - Checks all (including non finite-volume specific) geometry + - \par noTopology + 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 \ Specify an alternative mesh region. - - \par -writeSets \ - Reconstruct all cellSets and faceSets geometry and write to - postProcessing directory according to surfaceFormat - (e.g. vtk or ensight). Additionally reconstructs all pointSets and - writes as vtk format. + - \par -writeSurfaces + Reconstruct cellSets and faceSets of problem faces and write to + postProcessing directory. + + - \par -surfaceFormat + 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 used to write the pointSets + e.g. vtk or ensight. + + - \par -nonOrthThreshold + Threshold in degrees for reporting non-orthogonality errors, + default: 70" + + - \par -skewThreshold + 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; - const bool writeSets = args.optionReadIfPresent("writeSets", surfaceFormat); + word surfaceFormat(vtkSurfaceWriter::typeName); + 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 - << " representation" - << " of all faceSets and cellSets." << nl << endl; + Info<< "Reconstructing and writing the problem points in " + << setFormat << " format" + << nl << endl; } @@ -156,20 +222,23 @@ int main(int argc, char *argv[]) ); } - - autoPtr surfWriter; - autoPtr setWriter; - if (writeSets) + autoPtr surfaceWriter; + if (writeSurfaces) { - surfWriter = surfaceWriter::New + surfaceWriter = surfaceWriter::New ( surfaceFormat, mesh.time().writeFormat(), mesh.time().writeCompression() ); + } + + autoPtr 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); } diff --git a/src/dynamicMesh/polyTopoChange/attachPolyTopoChanger/attachPolyTopoChanger.C b/src/dynamicMesh/polyTopoChange/attachPolyTopoChanger/attachPolyTopoChanger.C index 55ea0c35cf..6b209d9e55 100644 --- a/src/dynamicMesh/polyTopoChange/attachPolyTopoChanger/attachPolyTopoChanger.C +++ b/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_); } diff --git a/src/functionObjects/utilities/Make/options b/src/functionObjects/utilities/Make/options index a90f0b5aef..c4e53ecdaf 100644 --- a/src/functionObjects/utilities/Make/options +++ b/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 diff --git a/src/functionObjects/utilities/checkMesh/checkMesh.C b/src/functionObjects/utilities/checkMesh/checkMesh.C index 83eed77a8d..54d26eb090 100644 --- a/src/functionObjects/utilities/checkMesh/checkMesh.C +++ b/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 + ( + "surfaceFormat", + vtkSurfaceWriter::typeName + ); + } + + writeSets_ = dict.lookupOrDefault("writeSets", false); + if (writeSets_) + { + setFormat_ = dict.lookupOrDefault + ( + "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 surfWriter; + if (writeSurfaces_) + { + surfWriter = surfaceWriter::New + ( + surfaceFormat_, + mesh_.time().writeFormat(), + mesh_.time().writeCompression() + ); + } + + autoPtr 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 { diff --git a/src/functionObjects/utilities/checkMesh/checkMesh.H b/src/functionObjects/utilities/checkMesh/checkMesh.H index 5dae115d37..da53736d93 100644 --- a/src/functionObjects/utilities/checkMesh/checkMesh.H +++ b/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 { diff --git a/src/meshCheck/checkGeometry.C b/src/meshCheck/checkGeometry.C index fd903cee67..3bcc0d8b5b 100644 --- a/src/meshCheck/checkGeometry.C +++ b/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& surfWriter, const autoPtr& 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++; diff --git a/src/meshCheck/meshCheck.H b/src/meshCheck/meshCheck.H index 566e7c3b32..aad6cae854 100644 --- a/src/meshCheck/meshCheck.H +++ b/src/meshCheck/meshCheck.H @@ -81,6 +81,8 @@ namespace meshCheck ( const polyMesh& mesh, const bool allGeometry, + const scalar nonOrthThreshold, + const scalar skewThreshold, const autoPtr&, const autoPtr& ); diff --git a/src/meshCheck/polyMeshCheck/polyMeshCheck.C b/src/meshCheck/polyMeshCheck/polyMeshCheck.C index f58943fe28..6bcbdba5bf 100644 --- a/src/meshCheck/polyMeshCheck/polyMeshCheck.C +++ b/src/meshCheck/polyMeshCheck/polyMeshCheck.C @@ -281,6 +281,7 @@ Foam::tmp 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) { diff --git a/src/meshCheck/polyMeshCheck/polyMeshCheck.H b/src/meshCheck/polyMeshCheck/polyMeshCheck.H index 82a2d60e7a..f4ac5dfa5a 100644 --- a/src/meshCheck/polyMeshCheck/polyMeshCheck.H +++ b/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 ); diff --git a/src/meshCheck/polyMeshCheck/polyMeshCheckQuality.C b/src/meshCheck/polyMeshCheck/polyMeshCheckQuality.C index 3a59f78c34..39e0c73bad 100644 --- a/src/meshCheck/polyMeshCheck/polyMeshCheckQuality.C +++ b/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(); diff --git a/src/meshCheck/primitiveMeshCheck/primitiveMeshCheck.C b/src/meshCheck/primitiveMeshCheck/primitiveMeshCheck.C index c04e09831b..b20a97bcc8 100644 --- a/src/meshCheck/primitiveMeshCheck/primitiveMeshCheck.C +++ b/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;