zhyang-dev/OpenFOAM-12
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

dynamicMeshCheck merged into meshCheck

Browse Source
This commit is contained in:
Henry Weller
2023-10-20 17:46:04 +01:00
parent a5f0b5310a
commit 6628666770
37 changed files with 6267 additions and 1510 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C
View File

@ -53,8 +53,7 @@ Description
#include "polyAddFace.H"
#include "combineFaces.H"
#include "removePoints.H"
#include "polyMeshCheck.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "polyTopoChangeMap.H"
#include "unitConversion.H"
#include "motionSmoother.H"

2
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
View File

@ -47,7 +47,7 @@ Description
#include "snapParameters.H"
#include "layerParameters.H"
#include "faceSet.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "polyTopoChange.H"
#include "cellModeller.H"
#include "uindirectPrimitivePatch.H"

6
applications/utilities/mesh/manipulation/checkMesh/checkMesh.C
View File

@ -61,9 +61,7 @@ Usage
#include "vtkSetWriter.H"
#include "IOdictionary.H"
#include "checkTools.H"
#include "checkTopology.H"
#include "checkGeometry.H"
#include "meshCheck.H"
#include "checkMeshQuality.H"
using namespace Foam;
@ -196,7 +194,7 @@ int main(int argc, char *argv[])
// Reconstruct globalMeshData
mesh.globalData();
meshTools::printMeshStats(mesh, allTopology);
meshCheck::printMeshStats(mesh, allTopology);
label nFailedChecks = 0;

5
applications/utilities/mesh/manipulation/checkMesh/checkMeshQuality.C
View File

@ -27,9 +27,8 @@ License
#include "polyMesh.H"
#include "cellSet.H"
#include "faceSet.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "surfaceWriter.H"
#include "checkTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -58,7 +57,7 @@ Foam::label Foam::checkMeshQuality
faces.write();
if (writer.valid())
{
meshTools::mergeAndWrite(writer(), faces);
meshCheck::mergeAndWrite(writer(), faces);
}
}
}

2
src/dynamicMesh/motionSmoother/badQualityToCell/badQualityToCell.C
View File

@ -25,7 +25,7 @@ License
#include "badQualityToCell.H"
#include "polyMesh.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "faceSet.H"
#include "addToRunTimeSelectionTable.H"

2
src/dynamicMesh/motionSmoother/badQualityToFace/badQualityToFace.C
View File

@ -25,7 +25,7 @@ License
#include "badQualityToFace.H"
#include "polyMesh.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "faceSet.H"
#include "addToRunTimeSelectionTable.H"

2
src/dynamicMesh/motionSmoother/motionSmootherAlgo.C
View File

@ -31,7 +31,7 @@ License
#include "pointConstraints.H"
#include "syncTools.H"
#include "meshTools.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

2
src/dynamicMesh/polyTopoChange/attachPolyTopoChanger/attachPolyTopoChanger.C
View File

@ -24,7 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "attachPolyTopoChanger.H"
#include "polyMeshCheck.H"
#include "meshCheck.H"
#include "polyTopoChange.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

2
src/dynamicMesh/polyTopoChange/polyTopoChange/edgeCollapser.C
View File

@ -31,7 +31,7 @@ License
#include "PointEdgeWave.H"
#include "globalIndex.H"
#include "removePoints.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "OFstream.H"
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //

2
src/functionObjects/utilities/checkMesh/checkMesh.C
View File

@ -25,7 +25,7 @@ License
#include "checkMesh.H"
#include "fvMesh.H"
#include "polyMeshCheck.H"
#include "meshCheck.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

2
src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C
View File

@ -29,7 +29,7 @@ License
#include "removePoints.H"
#include "faceSet.H"
#include "Time.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

2
src/mesh/snappyHexMesh/meshRefinement/meshRefinementProblemCells.C
View File

@ -33,7 +33,7 @@ License
#include "indirectPrimitivePatch.H"
#include "cellSet.H"
#include "searchableSurfaces.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "IOmanip.H"
#include "unitConversion.H"
#include "snappySnapDriver.H"

2
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
View File

@ -33,7 +33,7 @@ Description
#include "meshRefinement.H"
#include "removePoints.H"
#include "pointFields.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "unitConversion.H"
#include "pointSet.H"
#include "faceSet.H"

2
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
View File

@ -28,7 +28,7 @@ Description
#include "snappySnapDriver.H"
#include "motionSmoother.H"
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "polyTopoChange.H"
#include "syncTools.H"
#include "fvMesh.H"

13
src/meshCheck/Make/files
View File

@ -1,16 +1,13 @@
primitiveMeshCheck/primitiveMeshCheck.C
primitiveMeshCheck/primitiveMeshCheckPointNearness.C
primitiveMeshCheck/primitiveMeshCheckEdgeLength.C
primitiveMeshCheck/primitiveMeshTools.C
polyMeshCheck/polyMeshCheck.C
polyMeshCheck/polyMeshTools.C
polyMeshCheck/polyMeshCheckQuality.C
dynamicMeshCheck/dynamicMeshCheck.C
dynamicMeshCheck/motionSmootherAlgoCheck.C
checkTools/checkTools.C
checkTools/checkTopology.C
checkTools/checkGeometry.C
mergeAndWrite/mergeAndWrite.C
checkTopology.C
checkGeometry.C
checkMesh.C
LIB = $(FOAM_LIBBIN)/libmeshCheck

40
src/meshCheck/checkTools/checkGeometry.C → src/meshCheck/checkGeometry.C
View File

@ -23,8 +23,8 @@ License
\*---------------------------------------------------------------------------*/
#include "meshCheck.H"
#include "PatchTools.H"
#include "checkGeometry.H"
#include "polyMeshCheck.H"
#include "cellSet.H"
#include "faceSet.H"
@ -39,7 +39,7 @@ License
#include "writeFile.H"
#include "nonConformalCyclicPolyPatch.H"
#include "checkTools.H"
#include "mergeAndWrite.H"
#include "Time.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -608,7 +608,7 @@ Foam::label Foam::meshCheck::checkGeometry
nonAlignedPoints.write();
if (setWriter.valid())
{
meshTools::mergeAndWrite(setWriter, nonAlignedPoints);
meshCheck::mergeAndWrite(setWriter, nonAlignedPoints);
}
}
}
@ -643,7 +643,7 @@ Foam::label Foam::meshCheck::checkGeometry
cells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), cells);
meshCheck::mergeAndWrite(surfWriter(), cells);
}
}
}
@ -659,7 +659,7 @@ Foam::label Foam::meshCheck::checkGeometry
aspectCells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), aspectCells);
meshCheck::mergeAndWrite(surfWriter(), aspectCells);
}
}
}
@ -680,7 +680,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -702,7 +702,7 @@ Foam::label Foam::meshCheck::checkGeometry
cells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), cells);
meshCheck::mergeAndWrite(surfWriter(), cells);
}
}
}
@ -725,7 +725,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -747,7 +747,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -769,7 +769,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -793,7 +793,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -826,7 +826,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -851,7 +851,7 @@ Foam::label Foam::meshCheck::checkGeometry
points.write();
if (setWriter.valid())
{
meshTools::mergeAndWrite(setWriter, points);
meshCheck::mergeAndWrite(setWriter, points);
}
}
}
@ -874,7 +874,7 @@ Foam::label Foam::meshCheck::checkGeometry
nearPoints.write();
if (setWriter.valid())
{
meshTools::mergeAndWrite(setWriter, nearPoints);
meshCheck::mergeAndWrite(setWriter, nearPoints);
}
}
}
@ -898,7 +898,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -921,7 +921,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -942,7 +942,7 @@ Foam::label Foam::meshCheck::checkGeometry
cells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), cells);
meshCheck::mergeAndWrite(surfWriter(), cells);
}
}
}
@ -962,7 +962,7 @@ Foam::label Foam::meshCheck::checkGeometry
cells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), cells);
meshCheck::mergeAndWrite(surfWriter(), cells);
}
}
}
@ -983,7 +983,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -1004,7 +1004,7 @@ Foam::label Foam::meshCheck::checkGeometry
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}

4
src/meshCheck/dynamicMeshCheck/motionSmootherAlgoCheck.C → src/meshCheck/checkMesh.C
View File

@ -23,7 +23,8 @@ License
\*---------------------------------------------------------------------------*/
#include "dynamicMeshCheck.H"
#include "meshCheck.H"
#include "polyMeshCheck.H"
#include "IOmanip.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
@ -49,6 +50,7 @@ bool Foam::meshCheck::checkMesh
);
}
bool Foam::meshCheck::checkMesh
(
const bool report,

62
src/meshCheck/checkTools/checkTopology.H
View File

@ -1,62 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Tools for checking the mesh
SourceFiles
checkTopology.C
\*---------------------------------------------------------------------------*/
#ifndef checkTopology_H
#define checkTopology_H
#include "label.H"
#include "autoPtr.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
class polyMesh;
class surfaceWriter;
class setWriter;
namespace meshCheck
{
label checkTopology
(
const polyMesh& mesh,
const bool allTopology,
const autoPtr<surfaceWriter>& surfWriter,
const autoPtr<Foam::setWriter>& setWriter
);
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

22
src/meshCheck/checkTools/checkTopology.C → src/meshCheck/checkTopology.C
View File

@ -23,7 +23,7 @@ License
\*---------------------------------------------------------------------------*/
#include "checkTopology.H"
#include "meshCheck.H"
#include "polyMeshCheck.H"
#include "Time.H"
#include "regionSplit.H"
@ -34,7 +34,7 @@ License
#include "emptyPolyPatch.H"
#include "processorPolyPatch.H"
#include "surfaceWriter.H"
#include "checkTools.H"
#include "mergeAndWrite.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -131,7 +131,7 @@ Foam::label Foam::meshCheck::checkTopology
cells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), cells);
meshCheck::mergeAndWrite(surfWriter(), cells);
}
}
@ -156,7 +156,7 @@ Foam::label Foam::meshCheck::checkTopology
points.write();
if (setWriter.valid())
{
meshTools::mergeAndWrite(setWriter, points);
meshCheck::mergeAndWrite(setWriter, points);
}
}
}
@ -178,7 +178,7 @@ Foam::label Foam::meshCheck::checkTopology
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -198,7 +198,7 @@ Foam::label Foam::meshCheck::checkTopology
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -219,7 +219,7 @@ Foam::label Foam::meshCheck::checkTopology
cells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), cells);
meshCheck::mergeAndWrite(surfWriter(), cells);
}
}
@ -243,7 +243,7 @@ Foam::label Foam::meshCheck::checkTopology
faces.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), faces);
meshCheck::mergeAndWrite(surfWriter(), faces);
}
}
}
@ -298,7 +298,7 @@ Foam::label Foam::meshCheck::checkTopology
oneCells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), oneCells);
meshCheck::mergeAndWrite(surfWriter(), oneCells);
}
}
@ -314,7 +314,7 @@ Foam::label Foam::meshCheck::checkTopology
twoCells.write();
if (surfWriter.valid())
{
meshTools::mergeAndWrite(surfWriter(), twoCells);
meshCheck::mergeAndWrite(surfWriter(), twoCells);
}
}
}
@ -459,7 +459,7 @@ Foam::label Foam::meshCheck::checkTopology
points.write();
if (setWriter.valid())
{
meshTools::mergeAndWrite(setWriter, points);
meshCheck::mergeAndWrite(setWriter, points);
}
}
}

251
src/meshCheck/dynamicMeshCheck/dynamicMeshCheck.H
View File

@ -1,251 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2011-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
PolyMesh checking routines. Checks various criteria for a mesh and supplied
geometry, with the mesh only used for topology. Coupled patch aware (i.e.,
coupled faces are treated as internal).
SourceFiles
dynamicMeshCheck.C
\*---------------------------------------------------------------------------*/
#ifndef dynamicMeshCheck_H
#define dynamicMeshCheck_H
#include "pointFields.H"
#include "HashSet.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace meshCheck
{
//- Check face non-orthogonality
bool checkFaceDotProduct
(
const bool report,
const scalar orthWarn,
const polyMesh&,
const vectorField& cellCentres,
const vectorField& faceAreas,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
);
//- Check face pyramid volumes
bool checkFacePyramids
(
const bool report,
const scalar minPyrVol,
const polyMesh&,
const vectorField& cellCentres,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet*
);
//- Check face tetrahedra volumes
bool checkFaceTets
(
const bool report,
const scalar minPyrVol,
const polyMesh&,
const vectorField& cellCentres,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet*
);
//- Check face skewness
bool checkFaceSkewness
(
const bool report,
const scalar internalSkew,
const scalar boundarySkew,
const polyMesh& mesh,
const pointField& points,
const vectorField& cellCentres,
const vectorField& faceCentres,
const vectorField& faceAreas,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
);
//- Check interpolation weights (0.5 for regular mesh)
bool checkFaceWeights
(
const bool report,
const scalar warnWeight,
const polyMesh& mesh,
const vectorField& cellCentres,
const vectorField& faceCentres,
const vectorField& faceAreas,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
);
//- Cell volume ratio of neighbouring cells (1 for regular mesh)
bool checkVolRatio
(
const bool report,
const scalar warnRatio,
const polyMesh& mesh,
const scalarField& cellVolumes,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
);
//- Check convexity of angles in a face. See primitiveMesh for explanation.
bool checkFaceAngles
(
const bool report,
const scalar maxDeg,
const polyMesh& mesh,
const vectorField& faceAreas,
const pointField& p,
const labelList& checkFaces,
labelHashSet* setPtr
);
// Check the difference between normals of individual face-triangles (from
// face-centre decomposition) and the cell-cell centre vector
bool checkFaceTwist
(
const bool report,
const scalar minTwist,
const polyMesh&,
const vectorField& cellCentres,
const vectorField& faceAreas,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check consecutive face-triangle (from face-centre decomposition) normals
bool checkTriangleTwist
(
const bool report,
const scalar minTwist,
const polyMesh&,
const vectorField& faceAreas,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check for face areas v.s. sum of face-triangle (from face-centre
// decomposition) areas
bool checkFaceFlatness
(
const bool report,
const scalar minFlatness,
const polyMesh&,
const vectorField& faceAreas,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check for small faces
bool checkFaceArea
(
const bool report,
const scalar minArea,
const polyMesh&,
const vectorField& faceAreas,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check the area of internal faces v.s. boundary faces
bool checkCellDeterminant
(
const bool report,
const scalar minDet,
const polyMesh&,
const vectorField& faceAreas,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check mesh with mesh settings in dict. Collects incorrect faces
// in set. Returns true if one or more faces in error.
// Parallel ok.
bool checkMesh
(
const bool report,
const polyMesh& mesh,
const dictionary& dict,
labelHashSet& wrongFaces
);
//- Check (subset of mesh) with mesh settings in dict.
// Collects incorrect faces in set. Returns true if one
// or more faces in error. Parallel ok.
bool checkMesh
(
const bool report,
const polyMesh& mesh,
const dictionary& dict,
const labelList& checkFaces,
labelHashSet& wrongFaces
);
//- Check (subset of mesh including baffles) with mesh settings
// in dict. Collects incorrect faces in set. Returns true if one
// or more faces in error. Parallel ok.
bool checkMesh
(
const bool report,
const polyMesh& mesh,
const dictionary& dict,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet& wrongFaces
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace meshCheck
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

14
src/meshCheck/checkTools/checkTools.C → src/meshCheck/mergeAndWrite/mergeAndWrite.C
View File

@ -23,7 +23,7 @@ License
\*---------------------------------------------------------------------------*/
#include "checkTools.H"
#include "mergeAndWrite.H"
#include "polyMesh.H"
#include "globalMeshData.H"
#include "hexMatcher.H"
@ -46,7 +46,7 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
void Foam::meshTools::printMeshStats
void Foam::meshCheck::printMeshStats
(
const polyMesh& mesh,
const bool allTopology
@ -209,7 +209,7 @@ void Foam::meshTools::printMeshStats
}
void Foam::meshTools::mergeAndWrite
void Foam::meshCheck::mergeAndWrite
(
const polyMesh& mesh,
const surfaceWriter& writer,
@ -261,7 +261,7 @@ void Foam::meshTools::mergeAndWrite
}
Foam::fileName Foam::meshTools::checkMeshOutputDir(const polyMesh& mesh)
Foam::fileName Foam::meshCheck::checkMeshOutputDir(const polyMesh& mesh)
{
return
mesh.time().globalPath()
@ -276,7 +276,7 @@ Foam::fileName Foam::meshTools::checkMeshOutputDir(const polyMesh& mesh)
}
void Foam::meshTools::mergeAndWrite
void Foam::meshCheck::mergeAndWrite
(
const surfaceWriter& writer,
const faceSet& set
@ -294,7 +294,7 @@ void Foam::meshTools::mergeAndWrite
}
void Foam::meshTools::mergeAndWrite
void Foam::meshCheck::mergeAndWrite
(
const surfaceWriter& writer,
const cellSet& set
@ -381,7 +381,7 @@ void Foam::meshTools::mergeAndWrite
}
void Foam::meshTools::mergeAndWrite
void Foam::meshCheck::mergeAndWrite
(
const setWriter& writer,
const pointSet& set

8
src/meshCheck/checkTools/checkTools.H → src/meshCheck/mergeAndWrite/mergeAndWrite.H
View File

@ -25,12 +25,12 @@ Description
Tools for checking the mesh
SourceFiles
checkTools.C
mergeAndWrite.C
\*---------------------------------------------------------------------------*/
#ifndef checkTools_H
#define checkTools_H
#ifndef mergeAndWrite_H
#define mergeAndWrite_H
#include "indirectPrimitivePatch.H"
@ -45,7 +45,7 @@ namespace Foam
class faceSet;
class cellSet;
namespace meshTools
namespace meshCheck
{
//- Output directory for sets and surfaces
fileName checkMeshOutputDir(const polyMesh& mesh);

57
src/meshCheck/checkTools/checkGeometry.H → src/meshCheck/meshCheck.H
View File

@ -22,20 +22,24 @@ License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Description
Routines for checking mesh geometry
Functions for checking mesh topology and geometry
SourceFiles
checkTopology.C
checkGeometry.C
\*---------------------------------------------------------------------------*/
#ifndef checkGeometry_H
#define checkGeometry_H
#ifndef meshCheck_H
#define meshCheck_H
#include "label.H"
#include "HashSet.H"
#include "labelVector.H"
#include "setWriter.H"
#include "labelPair.H"
#include "primitiveMeshCheck.H"
#include "polyMeshCheck.H"
#include "mergeAndWrite.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -63,6 +67,15 @@ namespace meshCheck
//- Check 0th vertex on coupled faces
bool checkCoupledPoints(const polyMesh&, const bool report, labelHashSet*);
//- Check the topology
label checkTopology
(
const polyMesh& mesh,
const bool allTopology,
const autoPtr<surfaceWriter>& surfWriter,
const autoPtr<Foam::setWriter>& setWriter
);
//- Check the geometry
label checkGeometry
(
@ -71,6 +84,42 @@ namespace meshCheck
const autoPtr<surfaceWriter>&,
const autoPtr<setWriter>&
);
//- Check mesh with mesh settings in dict. Collects incorrect faces
// in set. Returns true if one or more faces in error.
// Parallel ok.
bool checkMesh
(
const bool report,
const polyMesh& mesh,
const dictionary& dict,
labelHashSet& wrongFaces
);
//- Check (subset of mesh) with mesh settings in dict.
// Collects incorrect faces in set. Returns true if one
// or more faces in error. Parallel ok.
bool checkMesh
(
const bool report,
const polyMesh& mesh,
const dictionary& dict,
const labelList& checkFaces,
labelHashSet& wrongFaces
);
//- Check (subset of mesh including baffles) with mesh settings
// in dict. Collects incorrect faces in set. Returns true if one
// or more faces in error. Parallel ok.
bool checkMesh
(
const bool report,
const polyMesh& mesh,
const dictionary& dict,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet& wrongFaces
);
}
}

260
src/meshCheck/polyMeshCheck/polyMeshCheck.C
View File

@ -23,13 +23,261 @@ License
\*---------------------------------------------------------------------------*/
#include "primitiveMeshCheck.H"
#include "polyMeshCheck.H"
#include "polyMeshTools.H"
#include "unitConversion.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::scalarField> Foam::meshCheck::faceOrthogonality
(
const polyMesh& mesh,
const vectorField& areas,
const vectorField& cc
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tortho(new scalarField(mesh.nFaces(), 1.0));
scalarField& ortho = tortho.ref();
// Internal faces
forAll(nei, facei)
{
ortho[facei] = meshCheck::faceOrthogonality
(
cc[own[facei]],
cc[nei[facei]],
areas[facei]
);
}
// Coupled faces
pointField neighbourCc;
syncTools::swapBoundaryCellPositions(mesh, cc, neighbourCc);
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start() + i;
label bFacei = facei - mesh.nInternalFaces();
ortho[facei] = meshCheck::faceOrthogonality
(
cc[own[facei]],
neighbourCc[bFacei],
areas[facei]
);
}
}
}
return tortho;
}
Foam::tmp<Foam::scalarField> Foam::meshCheck::faceSkewness
(
const polyMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
scalarField& skew = tskew.ref();
forAll(nei, facei)
{
skew[facei] = meshCheck::faceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]],
cellCtrs[nei[facei]]
);
}
// Boundary faces: consider them to have only skewness error.
// (i.e. treat as if mirror cell on other side)
pointField neighbourCc;
syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neighbourCc);
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start() + i;
label bFacei = facei - mesh.nInternalFaces();
skew[facei] = meshCheck::faceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]],
neighbourCc[bFacei]
);
}
}
else
{
forAll(pp, i)
{
label facei = pp.start() + i;
skew[facei] = meshCheck::boundaryFaceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]]
);
}
}
}
return tskew;
}
Foam::tmp<Foam::scalarField> Foam::meshCheck::faceWeights
(
const polyMesh& mesh,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tweight(new scalarField(mesh.nFaces(), 1.0));
scalarField& weight = tweight.ref();
// Internal faces
forAll(nei, facei)
{
const point& fc = fCtrs[facei];
const vector& fa = fAreas[facei];
scalar dOwn = mag(fa & (fc-cellCtrs[own[facei]]));
scalar dNei = mag(fa & (cellCtrs[nei[facei]]-fc));
weight[facei] = min(dNei,dOwn)/(dNei+dOwn+vSmall);
}
// Coupled faces
pointField neiCc;
syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neiCc);
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start() + i;
label bFacei = facei - mesh.nInternalFaces();
const point& fc = fCtrs[facei];
const vector& fa = fAreas[facei];
scalar dOwn = mag(fa & (fc-cellCtrs[own[facei]]));
scalar dNei = mag(fa & (neiCc[bFacei]-fc));
weight[facei] = min(dNei,dOwn)/(dNei+dOwn+vSmall);
}
}
}
return tweight;
}
Foam::tmp<Foam::scalarField> Foam::meshCheck::volRatio
(
const polyMesh& mesh,
const scalarField& vol
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tratio(new scalarField(mesh.nFaces(), 1.0));
scalarField& ratio = tratio.ref();
// Internal faces
forAll(nei, facei)
{
scalar volOwn = vol[own[facei]];
scalar volNei = vol[nei[facei]];
ratio[facei] = min(volOwn,volNei)/(max(volOwn, volNei)+vSmall);
}
// Coupled faces
scalarField neiVol;
syncTools::swapBoundaryCellList(mesh, vol, neiVol);
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start() + i;
label bFacei = facei - mesh.nInternalFaces();
scalar volOwn = vol[own[facei]];
scalar volNei = neiVol[bFacei];
ratio[facei] = min(volOwn,volNei)/(max(volOwn, volNei)+vSmall);
}
}
}
return tratio;
}
bool Foam::meshCheck::checkFaceOrthogonality
(
const polyMesh& mesh,
@ -47,7 +295,7 @@ bool Foam::meshCheck::checkFaceOrthogonality
// Calculate orthogonality for all internal and coupled boundary faces
// (1 for uncoupled boundary faces)
tmp<scalarField> tortho = meshTools::faceOrthogonality
tmp<scalarField> tortho = meshCheck::faceOrthogonality
(
mesh,
fAreas,
@ -173,7 +421,7 @@ bool Foam::meshCheck::checkFaceSkewness
// Warn if the skew correction vector is more than skewWarning times
// larger than the face area vector
tmp<scalarField> tskew = meshTools::faceSkewness
tmp<scalarField> tskew = meshCheck::faceSkewness
(
mesh,
points,
@ -389,7 +637,7 @@ bool Foam::meshCheck::checkCellDeterminant
InfoInFunction << "Checking for under-determined cells" << endl;
}
tmp<scalarField> tcellDeterminant = meshTools::cellDeterminant
tmp<scalarField> tcellDeterminant = meshCheck::cellDeterminant
(
mesh,
meshD,
@ -473,7 +721,7 @@ bool Foam::meshCheck::checkFaceWeight
const vectorField& fAreas = mesh.faceAreas();
const vectorField& cellCtrs = mesh.cellCentres();
tmp<scalarField> tfaceWght = meshTools::faceWeights
tmp<scalarField> tfaceWght = meshCheck::faceWeights
(
mesh,
fCtrs,
@ -568,7 +816,7 @@ bool Foam::meshCheck::checkVolRatio
const scalarField& cellVols = mesh.cellVolumes();
tmp<scalarField> tvolRatio = meshTools::volRatio(mesh, cellVols);
tmp<scalarField> tvolRatio = meshCheck::volRatio(mesh, cellVols);
scalarField& volRatio = tvolRatio.ref();

192
src/meshCheck/polyMeshCheck/polyMeshCheck.H
View File

@ -48,6 +48,41 @@ namespace Foam
namespace meshCheck
{
//- Generate orthogonality field. (1 for fully orthogonal, < 1 for
// non-orthogonal)
tmp<scalarField> faceOrthogonality
(
const polyMesh& mesh,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate skewness field
tmp<scalarField> faceSkewness
(
const polyMesh& mesh,
const pointField& points,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate interpolation factors field
tmp<scalarField> faceWeights
(
const polyMesh& mesh,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate volume ratio field
tmp<scalarField> volRatio
(
const polyMesh& mesh,
const scalarField& vol
);
//- Check non-orthogonality
bool checkFaceOrthogonality
(
@ -98,6 +133,163 @@ namespace meshCheck
labelHashSet* setPtr = nullptr
);
//- Check face non-orthogonality
bool checkFaceDotProduct
(
const bool report,
const scalar orthWarn,
const polyMesh&,
const vectorField& cellCentres,
const vectorField& faceAreas,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
);
//- Check face pyramid volumes
bool checkFacePyramids
(
const bool report,
const scalar minPyrVol,
const polyMesh&,
const vectorField& cellCentres,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet*
);
//- Check face tetrahedra volumes
bool checkFaceTets
(
const bool report,
const scalar minPyrVol,
const polyMesh&,
const vectorField& cellCentres,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet*
);
//- Check face skewness
bool checkFaceSkewness
(
const bool report,
const scalar internalSkew,
const scalar boundarySkew,
const polyMesh& mesh,
const pointField& points,
const vectorField& cellCentres,
const vectorField& faceCentres,
const vectorField& faceAreas,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
);
//- Check interpolation weights (0.5 for regular mesh)
bool checkFaceWeights
(
const bool report,
const scalar warnWeight,
const polyMesh& mesh,
const vectorField& cellCentres,
const vectorField& faceCentres,
const vectorField& faceAreas,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
);
//- Cell volume ratio of neighbouring cells (1 for regular mesh)
bool checkVolRatio
(
const bool report,
const scalar warnRatio,
const polyMesh& mesh,
const scalarField& cellVolumes,
const labelList& checkFaces,
const List<labelPair>& baffles,
labelHashSet* setPtr
);
//- Check convexity of angles in a face. See primitiveMesh for explanation.
bool checkFaceAngles
(
const bool report,
const scalar maxDeg,
const polyMesh& mesh,
const vectorField& faceAreas,
const pointField& p,
const labelList& checkFaces,
labelHashSet* setPtr
);
// Check the difference between normals of individual face-triangles (from
// face-centre decomposition) and the cell-cell centre vector
bool checkFaceTwist
(
const bool report,
const scalar minTwist,
const polyMesh&,
const vectorField& cellCentres,
const vectorField& faceAreas,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check consecutive face-triangle (from face-centre decomposition) normals
bool checkTriangleTwist
(
const bool report,
const scalar minTwist,
const polyMesh&,
const vectorField& faceAreas,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check for face areas v.s. sum of face-triangle (from face-centre
// decomposition) areas
bool checkFaceFlatness
(
const bool report,
const scalar minFlatness,
const polyMesh&,
const vectorField& faceAreas,
const vectorField& faceCentres,
const pointField& p,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check for small faces
bool checkFaceArea
(
const bool report,
const scalar minArea,
const polyMesh&,
const vectorField& faceAreas,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check the area of internal faces v.s. boundary faces
bool checkCellDeterminant
(
const bool report,
const scalar minDet,
const polyMesh&,
const vectorField& faceAreas,
const labelList& checkFaces,
labelHashSet* setPtr
);
//- Check mesh topology for correctness.
// Returns false for no error.

12
src/meshCheck/dynamicMeshCheck/dynamicMeshCheck.C → src/meshCheck/polyMeshCheck/polyMeshCheckQuality.C
View File

@ -23,13 +23,13 @@ License
\*---------------------------------------------------------------------------*/
#include "dynamicMeshCheck.H"
#include "primitiveMeshCheck.H"
#include "polyMeshCheck.H"
#include "polyMeshTetDecomposition.H"
#include "pyramidPointFaceRef.H"
#include "tetPointRef.H"
#include "syncTools.H"
#include "unitConversion.H"
#include "primitiveMeshTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -819,7 +819,7 @@ bool Foam::meshCheck::checkFaceSkewness
if (mesh.isInternalFace(facei))
{
scalar skewness = meshTools::faceSkewness
scalar skewness = meshCheck::faceSkewness
(
mesh,
points,
@ -854,7 +854,7 @@ bool Foam::meshCheck::checkFaceSkewness
}
else if (patches[patches.whichPatch(facei)].coupled())
{
scalar skewness = meshTools::faceSkewness
scalar skewness = meshCheck::faceSkewness
(
mesh,
points,
@ -889,7 +889,7 @@ bool Foam::meshCheck::checkFaceSkewness
}
else
{
scalar skewness = meshTools::boundaryFaceSkewness
scalar skewness = meshCheck::boundaryFaceSkewness
(
mesh,
points,
@ -932,7 +932,7 @@ bool Foam::meshCheck::checkFaceSkewness
const point& ownCc = cellCentres[own[face0]];
const point& neiCc = cellCentres[own[face1]];
scalar skewness = meshTools::faceSkewness
scalar skewness = meshCheck::faceSkewness
(
mesh,
points,

282
src/meshCheck/polyMeshCheck/polyMeshTools.C
View File

@ -1,282 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "polyMeshTools.H"
#include "syncTools.H"
#include "pyramidPointFaceRef.H"
#include "primitiveMeshTools.H"
#include "polyMeshTools.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::scalarField> Foam::meshTools::faceOrthogonality
(
const polyMesh& mesh,
const vectorField& areas,
const vectorField& cc
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tortho(new scalarField(mesh.nFaces(), 1.0));
scalarField& ortho = tortho.ref();
// Internal faces
forAll(nei, facei)
{
ortho[facei] = meshTools::faceOrthogonality
(
cc[own[facei]],
cc[nei[facei]],
areas[facei]
);
}
// Coupled faces
pointField neighbourCc;
syncTools::swapBoundaryCellPositions(mesh, cc, neighbourCc);
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start() + i;
label bFacei = facei - mesh.nInternalFaces();
ortho[facei] = meshTools::faceOrthogonality
(
cc[own[facei]],
neighbourCc[bFacei],
areas[facei]
);
}
}
}
return tortho;
}
Foam::tmp<Foam::scalarField> Foam::meshTools::faceSkewness
(
const polyMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
scalarField& skew = tskew.ref();
forAll(nei, facei)
{
skew[facei] = meshTools::faceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]],
cellCtrs[nei[facei]]
);
}
// Boundary faces: consider them to have only skewness error.
// (i.e. treat as if mirror cell on other side)
pointField neighbourCc;
syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neighbourCc);
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start() + i;
label bFacei = facei - mesh.nInternalFaces();
skew[facei] = meshTools::faceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]],
neighbourCc[bFacei]
);
}
}
else
{
forAll(pp, i)
{
label facei = pp.start() + i;
skew[facei] = meshTools::boundaryFaceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]]
);
}
}
}
return tskew;
}
Foam::tmp<Foam::scalarField> Foam::meshTools::faceWeights
(
const polyMesh& mesh,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tweight(new scalarField(mesh.nFaces(), 1.0));
scalarField& weight = tweight.ref();
// Internal faces
forAll(nei, facei)
{
const point& fc = fCtrs[facei];
const vector& fa = fAreas[facei];
scalar dOwn = mag(fa & (fc-cellCtrs[own[facei]]));
scalar dNei = mag(fa & (cellCtrs[nei[facei]]-fc));
weight[facei] = min(dNei,dOwn)/(dNei+dOwn+vSmall);
}
// Coupled faces
pointField neiCc;
syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neiCc);
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start() + i;
label bFacei = facei - mesh.nInternalFaces();
const point& fc = fCtrs[facei];
const vector& fa = fAreas[facei];
scalar dOwn = mag(fa & (fc-cellCtrs[own[facei]]));
scalar dNei = mag(fa & (neiCc[bFacei]-fc));
weight[facei] = min(dNei,dOwn)/(dNei+dOwn+vSmall);
}
}
}
return tweight;
}
Foam::tmp<Foam::scalarField> Foam::meshTools::volRatio
(
const polyMesh& mesh,
const scalarField& vol
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const polyBoundaryMesh& pbm = mesh.boundaryMesh();
tmp<scalarField> tratio(new scalarField(mesh.nFaces(), 1.0));
scalarField& ratio = tratio.ref();
// Internal faces
forAll(nei, facei)
{
scalar volOwn = vol[own[facei]];
scalar volNei = vol[nei[facei]];
ratio[facei] = min(volOwn,volNei)/(max(volOwn, volNei)+vSmall);
}
// Coupled faces
scalarField neiVol;
syncTools::swapBoundaryCellList(mesh, vol, neiVol);
forAll(pbm, patchi)
{
const polyPatch& pp = pbm[patchi];
if (pp.coupled())
{
forAll(pp, i)
{
label facei = pp.start() + i;
label bFacei = facei - mesh.nInternalFaces();
scalar volOwn = vol[own[facei]];
scalar volNei = neiVol[bFacei];
ratio[facei] = min(volOwn,volNei)/(max(volOwn, volNei)+vSmall);
}
}
}
return tratio;
}
// ************************************************************************* //

97
src/meshCheck/polyMeshCheck/polyMeshTools.H
View File

@ -1,97 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Namespace
Foam::meshTools
Description
Collection of functions operating on polyMesh (mainly checks) so
that need access to patch information.
SourceFiles
polyMeshTools.C
\*---------------------------------------------------------------------------*/
#ifndef polyMeshTools_H
#define polyMeshTools_H
#include "polyMesh.H"
#include "primitiveMeshTools.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Namespace meshTools Declaration
\*---------------------------------------------------------------------------*/
namespace meshTools
{
//- Generate orthogonality field. (1 for fully orthogonal, < 1 for
// non-orthogonal)
tmp<scalarField> faceOrthogonality
(
const polyMesh& mesh,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate skewness field
tmp<scalarField> faceSkewness
(
const polyMesh& mesh,
const pointField& points,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate interpolation factors field
tmp<scalarField> faceWeights
(
const polyMesh& mesh,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate volume ratio field
tmp<scalarField> volRatio
(
const polyMesh& mesh,
const scalarField& vol
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

529
src/meshCheck/primitiveMeshCheck/primitiveMeshCheck.C
View File

@ -25,11 +25,11 @@ License
#include "primitiveMeshCheck.H"
#include "pyramidPointFaceRef.H"
#include "PackedBoolList.H"
#include "ListOps.H"
#include "unitConversion.H"
#include "SortableList.H"
#include "EdgeMap.H"
#include "primitiveMeshTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -42,6 +42,525 @@ Foam::scalar Foam::meshCheck::planarCosAngle = 1.0e-6;
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::meshCheck::faceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label facei,
const point& ownCc,
const point& neiCc
)
{
vector Cpf = fCtrs[facei] - ownCc;
vector d = neiCc - ownCc;
// Skewness vector
vector sv =
Cpf
- ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + rootVSmall))*d;
vector svHat = sv/(mag(sv) + rootVSmall);
// Normalisation distance calculated as the approximate distance
// from the face centre to the edge of the face in the direction
// of the skewness
scalar fd = 0.2*mag(d) + rootVSmall;
const face& f = mesh.faces()[facei];
forAll(f, pi)
{
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
}
// Normalised skewness
return mag(sv)/fd;
}
Foam::scalar Foam::meshCheck::boundaryFaceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label facei,
const point& ownCc
)
{
vector Cpf = fCtrs[facei] - ownCc;
vector normal = fAreas[facei];
normal /= mag(normal) + rootVSmall;
vector d = normal*(normal & Cpf);
// Skewness vector
vector sv =
Cpf
- ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + rootVSmall))*d;
vector svHat = sv/(mag(sv) + rootVSmall);
// Normalisation distance calculated as the approximate distance
// from the face centre to the edge of the face in the direction
// of the skewness
scalar fd = 0.4*mag(d) + rootVSmall;
const face& f = mesh.faces()[facei];
forAll(f, pi)
{
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
}
// Normalised skewness
return mag(sv)/fd;
}
Foam::scalar Foam::meshCheck::faceOrthogonality
(
const point& ownCc,
const point& neiCc,
const vector& s
)
{
vector d = neiCc - ownCc;
return (d & s)/(mag(d)*mag(s) + rootVSmall);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::scalarField> Foam::meshCheck::faceOrthogonality
(
const primitiveMesh& mesh,
const vectorField& areas,
const vectorField& cc
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
tmp<scalarField> tortho(new scalarField(mesh.nInternalFaces()));
scalarField& ortho = tortho.ref();
// Internal faces
forAll(nei, facei)
{
ortho[facei] = faceOrthogonality
(
cc[own[facei]],
cc[nei[facei]],
areas[facei]
);
}
return tortho;
}
Foam::tmp<Foam::scalarField> Foam::meshCheck::faceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
scalarField& skew = tskew.ref();
forAll(nei, facei)
{
skew[facei] = faceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]],
cellCtrs[nei[facei]]
);
}
// Boundary faces: consider them to have only skewness error.
// (i.e. treat as if mirror cell on other side)
for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
{
skew[facei] = boundaryFaceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]]
);
}
return tskew;
}
void Foam::meshCheck::facePyramidVolume
(
const primitiveMesh& mesh,
const pointField& points,
const vectorField& ctrs,
scalarField& ownPyrVol,
scalarField& neiPyrVol
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const faceList& f = mesh.faces();
ownPyrVol.setSize(mesh.nFaces());
neiPyrVol.setSize(mesh.nInternalFaces());
forAll(f, facei)
{
// Create the owner pyramid
ownPyrVol[facei] = -pyramidPointFaceRef
(
f[facei],
ctrs[own[facei]]
).mag(points);
if (mesh.isInternalFace(facei))
{
// Create the neighbour pyramid - it will have positive volume
neiPyrVol[facei] = pyramidPointFaceRef
(
f[facei],
ctrs[nei[facei]]
).mag(points);
}
}
}
void Foam::meshCheck::cellClosedness
(
const primitiveMesh& mesh,
const Vector<label>& meshD,
const vectorField& areas,
const scalarField& vols,
scalarField& openness,
scalarField& aratio
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
// Loop through cell faces and sum up the face area vectors for each cell.
// This should be zero in all vector components
vectorField sumClosed(mesh.nCells(), Zero);
vectorField sumMagClosed(mesh.nCells(), Zero);
forAll(own, facei)
{
// Add to owner
sumClosed[own[facei]] += areas[facei];
sumMagClosed[own[facei]] += cmptMag(areas[facei]);
}
forAll(nei, facei)
{
// Subtract from neighbour
sumClosed[nei[facei]] -= areas[facei];
sumMagClosed[nei[facei]] += cmptMag(areas[facei]);
}
label nDims = 0;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
nDims++;
}
}
// Check the sums
openness.setSize(mesh.nCells());
aratio.setSize(mesh.nCells());
forAll(sumClosed, celli)
{
scalar maxOpenness = 0;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
maxOpenness = max
(
maxOpenness,
mag(sumClosed[celli][cmpt])
/(sumMagClosed[celli][cmpt] + rootVSmall)
);
}
openness[celli] = maxOpenness;
// Calculate the aspect ration as the maximum of Cartesian component
// aspect ratio to the total area hydraulic area aspect ratio
scalar minCmpt = vGreat;
scalar maxCmpt = -vGreat;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
minCmpt = min(minCmpt, sumMagClosed[celli][dir]);
maxCmpt = max(maxCmpt, sumMagClosed[celli][dir]);
}
}
scalar aspectRatio = maxCmpt/(minCmpt + rootVSmall);
if (nDims == 3)
{
scalar v = max(rootVSmall, vols[celli]);
aspectRatio = max
(
aspectRatio,
1.0/6.0*cmptSum(sumMagClosed[celli])/pow(v, 2.0/3.0)
);
}
aratio[celli] = aspectRatio;
}
}
Foam::tmp<Foam::scalarField> Foam::meshCheck::faceConcavity
(
const scalar maxSin,
const primitiveMesh& mesh,
const pointField& p,
const vectorField& faceAreas
)
{
const faceList& fcs = mesh.faces();
vectorField faceNormals(faceAreas);
faceNormals /= mag(faceNormals) + rootVSmall;
tmp<scalarField> tfaceAngles(new scalarField(mesh.nFaces()));
scalarField& faceAngles = tfaceAngles.ref();
forAll(fcs, facei)
{
const face& f = fcs[facei];
// Get edge from f[0] to f[size-1];
vector ePrev(p[f.first()] - p[f.last()]);
scalar magEPrev = mag(ePrev);
ePrev /= magEPrev + rootVSmall;
scalar maxEdgeSin = 0.0;
forAll(f, fp0)
{
// Get vertex after fp
label fp1 = f.fcIndex(fp0);
// Normalised vector between two consecutive points
vector e10(p[f[fp1]] - p[f[fp0]]);
scalar magE10 = mag(e10);
e10 /= magE10 + rootVSmall;
if (magEPrev > small && magE10 > small)
{
vector edgeNormal = ePrev ^ e10;
scalar magEdgeNormal = mag(edgeNormal);
if (magEdgeNormal < maxSin)
{
// Edges (almost) aligned -> face is ok.
}
else
{
// Check normal
edgeNormal /= magEdgeNormal;
if ((edgeNormal & faceNormals[facei]) < small)
{
maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
}
}
}
ePrev = e10;
magEPrev = magE10;
}
faceAngles[facei] = maxEdgeSin;
}
return tfaceAngles;
}
Foam::tmp<Foam::scalarField> Foam::meshCheck::faceFlatness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& faceAreas
)
{
const faceList& fcs = mesh.faces();
// Areas are calculated as the sum of areas. (see
// primitiveMeshFaceCentresAndAreas.C)
scalarField magAreas(mag(faceAreas));
tmp<scalarField> tfaceFlatness(new scalarField(mesh.nFaces(), 1.0));
scalarField& faceFlatness = tfaceFlatness.ref();
forAll(fcs, facei)
{
const face& f = fcs[facei];
if (f.size() > 3 && magAreas[facei] > rootVSmall)
{
const point& fc = fCtrs[facei];
// Calculate the sum of magnitude of areas and compare to magnitude
// of sum of areas.
scalar sumA = 0.0;
forAll(f, fp)
{
const point& thisPoint = p[f[fp]];
const point& nextPoint = p[f.nextLabel(fp)];
// Triangle around fc.
vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
sumA += mag(n);
}
faceFlatness[facei] = magAreas[facei]/(sumA + rootVSmall);
}
}
return tfaceFlatness;
}
Foam::tmp<Foam::scalarField> Foam::meshCheck::cellDeterminant
(
const primitiveMesh& mesh,
const Vector<label>& meshD,
const vectorField& faceAreas,
const PackedBoolList& internalOrCoupledFace
)
{
// Determine number of dimensions and (for 2D) missing dimension
label nDims = 0;
label twoD = -1;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
nDims++;
}
else
{
twoD = dir;
}
}
tmp<scalarField> tcellDeterminant(new scalarField(mesh.nCells()));
scalarField& cellDeterminant = tcellDeterminant.ref();
const cellList& c = mesh.cells();
if (nDims == 1)
{
cellDeterminant = 1.0;
}
else
{
forAll(c, celli)
{
const labelList& curFaces = c[celli];
// Calculate local normalisation factor
scalar avgArea = 0;
label nInternalFaces = 0;
forAll(curFaces, i)
{
if (internalOrCoupledFace[curFaces[i]])
{
avgArea += mag(faceAreas[curFaces[i]]);
nInternalFaces++;
}
}
if (nInternalFaces == 0)
{
cellDeterminant[celli] = 0;
}
else
{
avgArea /= nInternalFaces;
symmTensor areaTensor(Zero);
forAll(curFaces, i)
{
if (internalOrCoupledFace[curFaces[i]])
{
areaTensor += sqr(faceAreas[curFaces[i]]/avgArea);
}
}
if (nDims == 2)
{
// Add the missing eigenvector (such that it does not
// affect the determinant)
if (twoD == 0)
{
areaTensor.xx() = 1;
}
else if (twoD == 1)
{
areaTensor.yy() = 1;
}
else
{
areaTensor.zz() = 1;
}
}
cellDeterminant[celli] = mag(det(areaTensor));
}
}
}
return tcellDeterminant;
}
bool Foam::meshCheck::checkClosedBoundary
(
const primitiveMesh& mesh,
@ -149,7 +668,7 @@ bool Foam::meshCheck::checkClosedCells
scalarField openness;
scalarField aspectRatio;
meshTools::cellClosedness
meshCheck::cellClosedness
(
mesh,
meshD,
@ -374,7 +893,7 @@ bool Foam::meshCheck::checkFacePyramids
scalarField ownPyrVol;
scalarField neiPyrVol;
meshTools::facePyramidVolume
meshCheck::facePyramidVolume
(
mesh,
points,
@ -461,7 +980,7 @@ bool Foam::meshCheck::checkFaceAngles
const pointField& points = mesh.points();
const vectorField& faceAreas = mesh.faceAreas();
tmp<scalarField> tfaceAngles = meshTools::faceConcavity
tmp<scalarField> tfaceAngles = meshCheck::faceConcavity
(
maxSin,
mesh,
@ -542,7 +1061,7 @@ bool Foam::meshCheck::checkFaceFlatness
const vectorField& faceAreas = mesh.faceAreas();
const faceList& fcs = mesh.faces();
tmp<scalarField> tfaceFlatness = meshTools::faceFlatness
tmp<scalarField> tfaceFlatness = meshCheck::faceFlatness
(
mesh,
points,

114
src/meshCheck/primitiveMeshCheck/primitiveMeshCheck.H
View File

@ -130,6 +130,120 @@ namespace meshCheck
// Geometric checks
//- Generate non-orthogonality field (internal faces only)
tmp<scalarField> faceOrthogonality
(
const primitiveMesh& mesh,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate face pyramid volume fields
void facePyramidVolume
(
const primitiveMesh& mesh,
const pointField& points,
const vectorField& cellCtrs,
scalarField& ownPyrVol,
scalarField& neiPyrVol
);
//- Generate skewness field
tmp<scalarField> faceSkewness
(
const primitiveMesh& mesh,
const pointField& points,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate cell openness and cell aspect ratio field
void cellClosedness
(
const primitiveMesh& mesh,
const Vector<label>& meshD,
const vectorField& areas,
const scalarField& vols,
scalarField& openness,
scalarField& aratio
);
//- Generate face concavity field. Returns per face the (sin of the)
// most concave angle between two consecutive edges
tmp<scalarField> faceConcavity
(
const scalar maxSin,
const primitiveMesh& mesh,
const pointField& p,
const vectorField& faceAreas
);
//- Generate face flatness field. Compares the individual triangles'
// normals against the face average normal. Between 0 (fully warped)
// and 1 (fully flat)
tmp<scalarField> faceFlatness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& faceAreas
);
//- Generate edge alignment field. Is per face the minimum aligned edge
// (does not use edge addressing)
tmp<scalarField> edgeAlignment
(
const primitiveMesh& mesh,
const Vector<label>& directions,
const pointField& p
);
//- Generate cell determinant field
tmp<scalarField> cellDeterminant
(
const primitiveMesh& mesh,
const Vector<label>& directions,
const vectorField& faceAreas,
const PackedBoolList& internalOrCoupledFace
);
// Helpers: single face check
//- Skewness of single face
scalar faceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label facei,
const point& ownCc,
const point& neiCc
);
//- Skewness of single boundary face
scalar boundaryFaceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label facei,
const point& ownCc
);
//- Orthogonality of single face
scalar faceOrthogonality
(
const point& ownCc,
const point& neiCc,
const vector& s
);
//- Check boundary for closedness
bool checkClosedBoundary
(

551
src/meshCheck/primitiveMeshCheck/primitiveMeshTools.C
View File

@ -1,551 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "primitiveMeshTools.H"
#include "syncTools.H"
#include "pyramidPointFaceRef.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::scalar Foam::meshTools::faceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label facei,
const point& ownCc,
const point& neiCc
)
{
vector Cpf = fCtrs[facei] - ownCc;
vector d = neiCc - ownCc;
// Skewness vector
vector sv =
Cpf
- ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + rootVSmall))*d;
vector svHat = sv/(mag(sv) + rootVSmall);
// Normalisation distance calculated as the approximate distance
// from the face centre to the edge of the face in the direction
// of the skewness
scalar fd = 0.2*mag(d) + rootVSmall;
const face& f = mesh.faces()[facei];
forAll(f, pi)
{
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
}
// Normalised skewness
return mag(sv)/fd;
}
Foam::scalar Foam::meshTools::boundaryFaceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label facei,
const point& ownCc
)
{
vector Cpf = fCtrs[facei] - ownCc;
vector normal = fAreas[facei];
normal /= mag(normal) + rootVSmall;
vector d = normal*(normal & Cpf);
// Skewness vector
vector sv =
Cpf
- ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + rootVSmall))*d;
vector svHat = sv/(mag(sv) + rootVSmall);
// Normalisation distance calculated as the approximate distance
// from the face centre to the edge of the face in the direction
// of the skewness
scalar fd = 0.4*mag(d) + rootVSmall;
const face& f = mesh.faces()[facei];
forAll(f, pi)
{
fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei])));
}
// Normalised skewness
return mag(sv)/fd;
}
Foam::scalar Foam::meshTools::faceOrthogonality
(
const point& ownCc,
const point& neiCc,
const vector& s
)
{
vector d = neiCc - ownCc;
return (d & s)/(mag(d)*mag(s) + rootVSmall);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::tmp<Foam::scalarField> Foam::meshTools::faceOrthogonality
(
const primitiveMesh& mesh,
const vectorField& areas,
const vectorField& cc
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
tmp<scalarField> tortho(new scalarField(mesh.nInternalFaces()));
scalarField& ortho = tortho.ref();
// Internal faces
forAll(nei, facei)
{
ortho[facei] = faceOrthogonality
(
cc[own[facei]],
cc[nei[facei]],
areas[facei]
);
}
return tortho;
}
Foam::tmp<Foam::scalarField> Foam::meshTools::faceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
scalarField& skew = tskew.ref();
forAll(nei, facei)
{
skew[facei] = faceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]],
cellCtrs[nei[facei]]
);
}
// Boundary faces: consider them to have only skewness error.
// (i.e. treat as if mirror cell on other side)
for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
{
skew[facei] = boundaryFaceSkewness
(
mesh,
p,
fCtrs,
fAreas,
facei,
cellCtrs[own[facei]]
);
}
return tskew;
}
void Foam::meshTools::facePyramidVolume
(
const primitiveMesh& mesh,
const pointField& points,
const vectorField& ctrs,
scalarField& ownPyrVol,
scalarField& neiPyrVol
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
const faceList& f = mesh.faces();
ownPyrVol.setSize(mesh.nFaces());
neiPyrVol.setSize(mesh.nInternalFaces());
forAll(f, facei)
{
// Create the owner pyramid
ownPyrVol[facei] = -pyramidPointFaceRef
(
f[facei],
ctrs[own[facei]]
).mag(points);
if (mesh.isInternalFace(facei))
{
// Create the neighbour pyramid - it will have positive volume
neiPyrVol[facei] = pyramidPointFaceRef
(
f[facei],
ctrs[nei[facei]]
).mag(points);
}
}
}
void Foam::meshTools::cellClosedness
(
const primitiveMesh& mesh,
const Vector<label>& meshD,
const vectorField& areas,
const scalarField& vols,
scalarField& openness,
scalarField& aratio
)
{
const labelList& own = mesh.faceOwner();
const labelList& nei = mesh.faceNeighbour();
// Loop through cell faces and sum up the face area vectors for each cell.
// This should be zero in all vector components
vectorField sumClosed(mesh.nCells(), Zero);
vectorField sumMagClosed(mesh.nCells(), Zero);
forAll(own, facei)
{
// Add to owner
sumClosed[own[facei]] += areas[facei];
sumMagClosed[own[facei]] += cmptMag(areas[facei]);
}
forAll(nei, facei)
{
// Subtract from neighbour
sumClosed[nei[facei]] -= areas[facei];
sumMagClosed[nei[facei]] += cmptMag(areas[facei]);
}
label nDims = 0;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
nDims++;
}
}
// Check the sums
openness.setSize(mesh.nCells());
aratio.setSize(mesh.nCells());
forAll(sumClosed, celli)
{
scalar maxOpenness = 0;
for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
{
maxOpenness = max
(
maxOpenness,
mag(sumClosed[celli][cmpt])
/(sumMagClosed[celli][cmpt] + rootVSmall)
);
}
openness[celli] = maxOpenness;
// Calculate the aspect ration as the maximum of Cartesian component
// aspect ratio to the total area hydraulic area aspect ratio
scalar minCmpt = vGreat;
scalar maxCmpt = -vGreat;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
minCmpt = min(minCmpt, sumMagClosed[celli][dir]);
maxCmpt = max(maxCmpt, sumMagClosed[celli][dir]);
}
}
scalar aspectRatio = maxCmpt/(minCmpt + rootVSmall);
if (nDims == 3)
{
scalar v = max(rootVSmall, vols[celli]);
aspectRatio = max
(
aspectRatio,
1.0/6.0*cmptSum(sumMagClosed[celli])/pow(v, 2.0/3.0)
);
}
aratio[celli] = aspectRatio;
}
}
Foam::tmp<Foam::scalarField> Foam::meshTools::faceConcavity
(
const scalar maxSin,
const primitiveMesh& mesh,
const pointField& p,
const vectorField& faceAreas
)
{
const faceList& fcs = mesh.faces();
vectorField faceNormals(faceAreas);
faceNormals /= mag(faceNormals) + rootVSmall;
tmp<scalarField> tfaceAngles(new scalarField(mesh.nFaces()));
scalarField& faceAngles = tfaceAngles.ref();
forAll(fcs, facei)
{
const face& f = fcs[facei];
// Get edge from f[0] to f[size-1];
vector ePrev(p[f.first()] - p[f.last()]);
scalar magEPrev = mag(ePrev);
ePrev /= magEPrev + rootVSmall;
scalar maxEdgeSin = 0.0;
forAll(f, fp0)
{
// Get vertex after fp
label fp1 = f.fcIndex(fp0);
// Normalised vector between two consecutive points
vector e10(p[f[fp1]] - p[f[fp0]]);
scalar magE10 = mag(e10);
e10 /= magE10 + rootVSmall;
if (magEPrev > small && magE10 > small)
{
vector edgeNormal = ePrev ^ e10;
scalar magEdgeNormal = mag(edgeNormal);
if (magEdgeNormal < maxSin)
{
// Edges (almost) aligned -> face is ok.
}
else
{
// Check normal
edgeNormal /= magEdgeNormal;
if ((edgeNormal & faceNormals[facei]) < small)
{
maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
}
}
}
ePrev = e10;
magEPrev = magE10;
}
faceAngles[facei] = maxEdgeSin;
}
return tfaceAngles;
}
Foam::tmp<Foam::scalarField> Foam::meshTools::faceFlatness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& faceAreas
)
{
const faceList& fcs = mesh.faces();
// Areas are calculated as the sum of areas. (see
// primitiveMeshFaceCentresAndAreas.C)
scalarField magAreas(mag(faceAreas));
tmp<scalarField> tfaceFlatness(new scalarField(mesh.nFaces(), 1.0));
scalarField& faceFlatness = tfaceFlatness.ref();
forAll(fcs, facei)
{
const face& f = fcs[facei];
if (f.size() > 3 && magAreas[facei] > rootVSmall)
{
const point& fc = fCtrs[facei];
// Calculate the sum of magnitude of areas and compare to magnitude
// of sum of areas.
scalar sumA = 0.0;
forAll(f, fp)
{
const point& thisPoint = p[f[fp]];
const point& nextPoint = p[f.nextLabel(fp)];
// Triangle around fc.
vector n = 0.5*((nextPoint - thisPoint)^(fc - thisPoint));
sumA += mag(n);
}
faceFlatness[facei] = magAreas[facei]/(sumA + rootVSmall);
}
}
return tfaceFlatness;
}
Foam::tmp<Foam::scalarField> Foam::meshTools::cellDeterminant
(
const primitiveMesh& mesh,
const Vector<label>& meshD,
const vectorField& faceAreas,
const PackedBoolList& internalOrCoupledFace
)
{
// Determine number of dimensions and (for 2D) missing dimension
label nDims = 0;
label twoD = -1;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (meshD[dir] == 1)
{
nDims++;
}
else
{
twoD = dir;
}
}
tmp<scalarField> tcellDeterminant(new scalarField(mesh.nCells()));
scalarField& cellDeterminant = tcellDeterminant.ref();
const cellList& c = mesh.cells();
if (nDims == 1)
{
cellDeterminant = 1.0;
}
else
{
forAll(c, celli)
{
const labelList& curFaces = c[celli];
// Calculate local normalisation factor
scalar avgArea = 0;
label nInternalFaces = 0;
forAll(curFaces, i)
{
if (internalOrCoupledFace[curFaces[i]])
{
avgArea += mag(faceAreas[curFaces[i]]);
nInternalFaces++;
}
}
if (nInternalFaces == 0)
{
cellDeterminant[celli] = 0;
}
else
{
avgArea /= nInternalFaces;
symmTensor areaTensor(Zero);
forAll(curFaces, i)
{
if (internalOrCoupledFace[curFaces[i]])
{
areaTensor += sqr(faceAreas[curFaces[i]]/avgArea);
}
}
if (nDims == 2)
{
// Add the missing eigenvector (such that it does not
// affect the determinant)
if (twoD == 0)
{
areaTensor.xx() = 1;
}
else if (twoD == 1)
{
areaTensor.yy() = 1;
}
else
{
areaTensor.zz() = 1;
}
}
cellDeterminant[celli] = mag(det(areaTensor));
}
}
}
return tcellDeterminant;
}
// ************************************************************************* //

175
src/meshCheck/primitiveMeshCheck/primitiveMeshTools.H
View File

@ -1,175 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Copyright (C) 2012-2023 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Namespace
Foam::meshTools
Description
Collection of functions operating on primitiveMesh (mainly checks).
SourceFiles
primitiveMeshTools.C
\*---------------------------------------------------------------------------*/
#ifndef primitiveMeshTools_H
#define primitiveMeshTools_H
#include "primitiveMesh.H"
#include "PackedBoolList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Namespace meshTools Declaration
\*---------------------------------------------------------------------------*/
namespace meshTools
{
//- Generate non-orthogonality field (internal faces only)
tmp<scalarField> faceOrthogonality
(
const primitiveMesh& mesh,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate face pyramid volume fields
void facePyramidVolume
(
const primitiveMesh& mesh,
const pointField& points,
const vectorField& cellCtrs,
scalarField& ownPyrVol,
scalarField& neiPyrVol
);
//- Generate skewness field
tmp<scalarField> faceSkewness
(
const primitiveMesh& mesh,
const pointField& points,
const vectorField& fCtrs,
const vectorField& fAreas,
const vectorField& cellCtrs
);
//- Generate cell openness and cell aspect ratio field
void cellClosedness
(
const primitiveMesh& mesh,
const Vector<label>& meshD,
const vectorField& areas,
const scalarField& vols,
scalarField& openness,
scalarField& aratio
);
//- Generate face concavity field. Returns per face the (sin of the)
// most concave angle between two consecutive edges
tmp<scalarField> faceConcavity
(
const scalar maxSin,
const primitiveMesh& mesh,
const pointField& p,
const vectorField& faceAreas
);
//- Generate face flatness field. Compares the individual triangles'
// normals against the face average normal. Between 0 (fully warped)
// and 1 (fully flat)
tmp<scalarField> faceFlatness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& faceAreas
);
//- Generate edge alignment field. Is per face the minimum aligned edge
// (does not use edge addressing)
tmp<scalarField> edgeAlignment
(
const primitiveMesh& mesh,
const Vector<label>& directions,
const pointField& p
);
//- Generate cell determinant field
tmp<scalarField> cellDeterminant
(
const primitiveMesh& mesh,
const Vector<label>& directions,
const vectorField& faceAreas,
const PackedBoolList& internalOrCoupledFace
);
// Helpers: single face check
//- Skewness of single face
scalar faceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label facei,
const point& ownCc,
const point& neiCc
);
//- Skewness of single boundary face
scalar boundaryFaceSkewness
(
const primitiveMesh& mesh,
const pointField& p,
const vectorField& fCtrs,
const vectorField& fAreas,
const label facei,
const point& ownCc
);
//- Orthogonality of single face
scalar faceOrthogonality
(
const point& ownCc,
const point& neiCc,
const vector& s
);
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

42
tutorials/incompressibleFluid/cavity/constant/polyMesh/boundary Normal file
View File

@ -0,0 +1,42 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: dev
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class polyBoundaryMesh;
location "constant/polyMesh";
object boundary;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
3
(
movingWall
{
type wall;
inGroups List<word> 1(wall);
nFaces 20;
startFace 760;
}
fixedWalls
{
type wall;
inGroups List<word> 1(wall);
nFaces 60;
startFace 780;
}
frontAndBack
{
type empty;
inGroups List<word> 1(empty);
nFaces 800;
startFace 840;
}
)
// ************************************************************************* //

1663
tutorials/incompressibleFluid/cavity/constant/polyMesh/faces Normal file
View File

File diff suppressed because it is too large Load Diff

784
tutorials/incompressibleFluid/cavity/constant/polyMesh/neighbour Normal file
View File

@ -0,0 +1,784 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: dev
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class labelList;
note "nPoints: 882 nCells: 400 nFaces: 1640 nInternalFaces: 760";
location "constant/polyMesh";
object neighbour;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
760
(
1
20
2
21
3
22
4
23
5
24
6
25
7
26
8
27
9
28
10
29
11
30
12
31
13
32
14
33
15
34
16
35
17
36
18
37
19
38
39
21
40
22
41
23
42
24
43
25
44
26
45
27
46
28
47
29
48
30
49
31
50
32
51
33
52
34
53
35
54
36
55
37
56
38
57
39
58
59
41
60
42
61
43
62
44
63
45
64
46
65
47
66
48
67
49
68
50
69
51
70
52
71
53
72
54
73
55
74
56
75
57
76
58
77
59
78
79
61
80
62
81
63
82
64
83
65
84
66
85
67
86
68
87
69
88
70
89
71
90
72
91
73
92
74
93
75
94
76
95
77
96
78
97
79
98
99
81
100
82
101
83
102
84
103
85
104
86
105
87
106
88
107
89
108
90
109
91
110
92
111
93
112
94
113
95
114
96
115
97
116
98
117
99
118
119
101
120
102
121
103
122
104
123
105
124
106
125
107
126
108
127
109
128
110
129
111
130
112
131
113
132
114
133
115
134
116
135
117
136
118
137
119
138
139
121
140
122
141
123
142
124
143
125
144
126
145
127
146
128
147
129
148
130
149
131
150
132
151
133
152
134
153
135
154
136
155
137
156
138
157
139
158
159
141
160
142
161
143
162
144
163
145
164
146
165
147
166
148
167
149
168
150
169
151
170
152
171
153
172
154
173
155
174
156
175
157
176
158
177
159
178
179
161
180
162
181
163
182
164
183
165
184
166
185
167
186
168
187
169
188
170
189
171
190
172
191
173
192
174
193
175
194
176
195
177
196
178
197
179
198
199
181
200
182
201
183
202
184
203
185
204
186
205
187
206
188
207
189
208
190
209
191
210
192
211
193
212
194
213
195
214
196
215
197
216
198
217
199
218
219
201
220
202
221
203
222
204
223
205
224
206
225
207
226
208
227
209
228
210
229
211
230
212
231
213
232
214
233
215
234
216
235
217
236
218
237
219
238
239
221
240
222
241
223
242
224
243
225
244
226
245
227
246
228
247
229
248
230
249
231
250
232
251
233
252
234
253
235
254
236
255
237
256
238
257
239
258
259
241
260
242
261
243
262
244
263
245
264
246
265
247
266
248
267
249
268
250
269
251
270
252
271
253
272
254
273
255
274
256
275
257
276
258
277
259
278
279
261
280
262
281
263
282
264
283
265
284
266
285
267
286
268
287
269
288
270
289
271
290
272
291
273
292
274
293
275
294
276
295
277
296
278
297
279
298
299
281
300
282
301
283
302
284
303
285
304
286
305
287
306
288
307
289
308
290
309
291
310
292
311
293
312
294
313
295
314
296
315
297
316
298
317
299
318
319
301
320
302
321
303
322
304
323
305
324
306
325
307
326
308
327
309
328
310
329
311
330
312
331
313
332
314
333
315
334
316
335
317
336
318
337
319
338
339
321
340
322
341
323
342
324
343
325
344
326
345
327
346
328
347
329
348
330
349
331
350
332
351
333
352
334
353
335
354
336
355
337
356
338
357
339
358
359
341
360
342
361
343
362
344
363
345
364
346
365
347
366
348
367
349
368
350
369
351
370
352
371
353
372
354
373
355
374
356
375
357
376
358
377
359
378
379
361
380
362
381
363
382
364
383
365
384
366
385
367
386
368
387
369
388
370
389
371
390
372
391
373
392
374
393
375
394
376
395
377
396
378
397
379
398
399
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
)
// ************************************************************************* //

1664
tutorials/incompressibleFluid/cavity/constant/polyMesh/owner Normal file
View File

File diff suppressed because it is too large Load Diff

905
tutorials/incompressibleFluid/cavity/constant/polyMesh/points Normal file
View File

@ -0,0 +1,905 @@
/*--------------------------------*- C++ -*----------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
\\ / A nd | Version: dev
\\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
format ascii;
class vectorField;
location "constant/polyMesh";
object points;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
882
(
(0 0 0)
(0.005 0 0)
(0.01 0 0)
(0.015 0 0)
(0.02 0 0)
(0.025 0 0)
(0.03 0 0)
(0.035 0 0)
(0.04 0 0)
(0.045 0 0)
(0.05 0 0)
(0.055 0 0)
(0.06 0 0)
(0.065 0 0)
(0.07 0 0)
(0.075 0 0)
(0.08 0 0)
(0.085 0 0)
(0.09 0 0)
(0.095 0 0)
(0.1 0 0)
(0 0.005 0)
(0.005 0.005 0)
(0.01 0.005 0)
(0.015 0.005 0)
(0.02 0.005 0)
(0.025 0.005 0)
(0.03 0.005 0)
(0.035 0.005 0)
(0.04 0.005 0)
(0.045 0.005 0)
(0.05 0.005 0)
(0.055 0.005 0)
(0.06 0.005 0)
(0.065 0.005 0)
(0.07 0.005 0)
(0.075 0.005 0)
(0.08 0.005 0)
(0.085 0.005 0)
(0.09 0.005 0)
(0.095 0.005 0)
(0.1 0.005 0)
(0 0.01 0)
(0.005 0.01 0)
(0.01 0.01 0)
(0.015 0.01 0)
(0.02 0.01 0)
(0.025 0.01 0)
(0.03 0.01 0)
(0.035 0.01 0)
(0.04 0.01 0)
(0.045 0.01 0)
(0.05 0.01 0)
(0.055 0.01 0)
(0.06 0.01 0)
(0.065 0.01 0)
(0.07 0.01 0)
(0.075 0.01 0)
(0.08 0.01 0)
(0.085 0.01 0)
(0.09 0.01 0)
(0.095 0.01 0)
(0.1 0.01 0)
(0 0.015 0)
(0.005 0.015 0)
(0.01 0.015 0)
(0.015 0.015 0)
(0.02 0.015 0)
(0.025 0.015 0)
(0.03 0.015 0)
(0.035 0.015 0)
(0.04 0.015 0)
(0.045 0.015 0)
(0.05 0.015 0)
(0.055 0.015 0)
(0.06 0.015 0)
(0.065 0.015 0)
(0.07 0.015 0)
(0.075 0.015 0)
(0.08 0.015 0)
(0.085 0.015 0)
(0.09 0.015 0)
(0.095 0.015 0)
(0.1 0.015 0)
(0 0.02 0)
(0.005 0.02 0)
(0.01 0.02 0)
(0.015 0.02 0)
(0.02 0.02 0)
(0.025 0.02 0)
(0.03 0.02 0)
(0.035 0.02 0)
(0.04 0.02 0)
(0.045 0.02 0)
(0.05 0.02 0)
(0.055 0.02 0)
(0.06 0.02 0)
(0.065 0.02 0)
(0.07 0.02 0)
(0.075 0.02 0)
(0.08 0.02 0)
(0.085 0.02 0)
(0.09 0.02 0)
(0.095 0.02 0)
(0.1 0.02 0)
(0 0.025 0)
(0.005 0.025 0)
(0.01 0.025 0)
(0.015 0.025 0)
(0.02 0.025 0)
(0.025 0.025 0)
(0.03 0.025 0)
(0.035 0.025 0)
(0.04 0.025 0)
(0.045 0.025 0)
(0.05 0.025 0)
(0.055 0.025 0)
(0.06 0.025 0)
(0.065 0.025 0)
(0.07 0.025 0)
(0.075 0.025 0)
(0.08 0.025 0)
(0.085 0.025 0)
(0.09 0.025 0)
(0.095 0.025 0)
(0.1 0.025 0)
(0 0.03 0)
(0.005 0.03 0)
(0.01 0.03 0)
(0.015 0.03 0)
(0.02 0.03 0)
(0.025 0.03 0)
(0.03 0.03 0)
(0.035 0.03 0)
(0.04 0.03 0)
(0.045 0.03 0)
(0.05 0.03 0)
(0.055 0.03 0)
(0.06 0.03 0)
(0.065 0.03 0)
(0.07 0.03 0)
(0.075 0.03 0)
(0.08 0.03 0)
(0.085 0.03 0)
(0.09 0.03 0)
(0.095 0.03 0)
(0.1 0.03 0)
(0 0.035 0)
(0.005 0.035 0)
(0.01 0.035 0)
(0.015 0.035 0)
(0.02 0.035 0)
(0.025 0.035 0)
(0.03 0.035 0)
(0.035 0.035 0)
(0.04 0.035 0)
(0.045 0.035 0)
(0.05 0.035 0)
(0.055 0.035 0)
(0.06 0.035 0)
(0.065 0.035 0)
(0.07 0.035 0)
(0.075 0.035 0)
(0.08 0.035 0)
(0.085 0.035 0)
(0.09 0.035 0)
(0.095 0.035 0)
(0.1 0.035 0)
(0 0.04 0)
(0.005 0.04 0)
(0.01 0.04 0)
(0.015 0.04 0)
(0.02 0.04 0)
(0.025 0.04 0)
(0.03 0.04 0)
(0.035 0.04 0)
(0.04 0.04 0)
(0.045 0.04 0)
(0.05 0.04 0)
(0.055 0.04 0)
(0.06 0.04 0)
(0.065 0.04 0)
(0.07 0.04 0)
(0.075 0.04 0)
(0.08 0.04 0)
(0.085 0.04 0)
(0.09 0.04 0)
(0.095 0.04 0)
(0.1 0.04 0)
(0 0.045 0)
(0.005 0.045 0)
(0.01 0.045 0)
(0.015 0.045 0)
(0.02 0.045 0)
(0.025 0.045 0)
(0.03 0.045 0)
(0.035 0.045 0)
(0.04 0.045 0)
(0.045 0.045 0)
(0.05 0.045 0)
(0.055 0.045 0)
(0.06 0.045 0)
(0.065 0.045 0)
(0.07 0.045 0)
(0.075 0.045 0)
(0.08 0.045 0)
(0.085 0.045 0)
(0.09 0.045 0)
(0.095 0.045 0)
(0.1 0.045 0)
(0 0.05 0)
(0.005 0.05 0)
(0.01 0.05 0)
(0.015 0.05 0)
(0.02 0.05 0)
(0.025 0.05 0)
(0.03 0.05 0)
(0.035 0.05 0)
(0.04 0.05 0)
(0.045 0.05 0)
(0.05 0.05 0)
(0.055 0.05 0)
(0.06 0.05 0)
(0.065 0.05 0)
(0.07 0.05 0)
(0.075 0.05 0)
(0.08 0.05 0)
(0.085 0.05 0)
(0.09 0.05 0)
(0.095 0.05 0)
(0.1 0.05 0)
(0 0.055 0)
(0.005 0.055 0)
(0.01 0.055 0)
(0.015 0.055 0)
(0.02 0.055 0)
(0.025 0.055 0)
(0.03 0.055 0)
(0.035 0.055 0)
(0.04 0.055 0)
(0.045 0.055 0)
(0.05 0.055 0)
(0.055 0.055 0)
(0.06 0.055 0)
(0.065 0.055 0)
(0.07 0.055 0)
(0.075 0.055 0)
(0.08 0.055 0)
(0.085 0.055 0)
(0.09 0.055 0)
(0.095 0.055 0)
(0.1 0.055 0)
(0 0.06 0)
(0.005 0.06 0)
(0.01 0.06 0)
(0.015 0.06 0)
(0.02 0.06 0)
(0.025 0.06 0)
(0.03 0.06 0)
(0.035 0.06 0)
(0.04 0.06 0)
(0.045 0.06 0)
(0.05 0.06 0)
(0.055 0.06 0)
(0.06 0.06 0)
(0.065 0.06 0)
(0.07 0.06 0)
(0.075 0.06 0)
(0.08 0.06 0)
(0.085 0.06 0)
(0.09 0.06 0)
(0.095 0.06 0)
(0.1 0.06 0)
(0 0.065 0)
(0.005 0.065 0)
(0.01 0.065 0)
(0.015 0.065 0)
(0.02 0.065 0)
(0.025 0.065 0)
(0.03 0.065 0)
(0.035 0.065 0)
(0.04 0.065 0)
(0.045 0.065 0)
(0.05 0.065 0)
(0.055 0.065 0)
(0.06 0.065 0)
(0.065 0.065 0)
(0.07 0.065 0)
(0.075 0.065 0)
(0.08 0.065 0)
(0.085 0.065 0)
(0.09 0.065 0)
(0.095 0.065 0)
(0.1 0.065 0)
(0 0.07 0)
(0.005 0.07 0)
(0.01 0.07 0)
(0.015 0.07 0)
(0.02 0.07 0)
(0.025 0.07 0)
(0.03 0.07 0)
(0.035 0.07 0)
(0.04 0.07 0)
(0.045 0.07 0)
(0.05 0.07 0)
(0.055 0.07 0)
(0.06 0.07 0)
(0.065 0.07 0)
(0.07 0.07 0)
(0.075 0.07 0)
(0.08 0.07 0)
(0.085 0.07 0)
(0.09 0.07 0)
(0.095 0.07 0)
(0.1 0.07 0)
(0 0.075 0)
(0.005 0.075 0)
(0.01 0.075 0)
(0.015 0.075 0)
(0.02 0.075 0)
(0.025 0.075 0)
(0.03 0.075 0)
(0.035 0.075 0)
(0.04 0.075 0)
(0.045 0.075 0)
(0.05 0.075 0)
(0.055 0.075 0)
(0.06 0.075 0)
(0.065 0.075 0)
(0.07 0.075 0)
(0.075 0.075 0)
(0.08 0.075 0)
(0.085 0.075 0)
(0.09 0.075 0)
(0.095 0.075 0)
(0.1 0.075 0)
(0 0.08 0)
(0.005 0.08 0)
(0.01 0.08 0)
(0.015 0.08 0)
(0.02 0.08 0)
(0.025 0.08 0)
(0.03 0.08 0)
(0.035 0.08 0)
(0.04 0.08 0)
(0.045 0.08 0)
(0.05 0.08 0)
(0.055 0.08 0)
(0.06 0.08 0)
(0.065 0.08 0)
(0.07 0.08 0)
(0.075 0.08 0)
(0.08 0.08 0)
(0.085 0.08 0)
(0.09 0.08 0)
(0.095 0.08 0)
(0.1 0.08 0)
(0 0.085 0)
(0.005 0.085 0)
(0.01 0.085 0)
(0.015 0.085 0)
(0.02 0.085 0)
(0.025 0.085 0)
(0.03 0.085 0)
(0.035 0.085 0)
(0.04 0.085 0)
(0.045 0.085 0)
(0.05 0.085 0)
(0.055 0.085 0)
(0.06 0.085 0)
(0.065 0.085 0)
(0.07 0.085 0)
(0.075 0.085 0)
(0.08 0.085 0)
(0.085 0.085 0)
(0.09 0.085 0)
(0.095 0.085 0)
(0.1 0.085 0)
(0 0.09 0)
(0.005 0.09 0)
(0.01 0.09 0)
(0.015 0.09 0)
(0.02 0.09 0)
(0.025 0.09 0)
(0.03 0.09 0)
(0.035 0.09 0)
(0.04 0.09 0)
(0.045 0.09 0)
(0.05 0.09 0)
(0.055 0.09 0)
(0.06 0.09 0)
(0.065 0.09 0)
(0.07 0.09 0)
(0.075 0.09 0)
(0.08 0.09 0)
(0.085 0.09 0)
(0.09 0.09 0)
(0.095 0.09 0)
(0.1 0.09 0)
(0 0.095 0)
(0.005 0.095 0)
(0.01 0.095 0)
(0.015 0.095 0)
(0.02 0.095 0)
(0.025 0.095 0)
(0.03 0.095 0)
(0.035 0.095 0)
(0.04 0.095 0)
(0.045 0.095 0)
(0.05 0.095 0)
(0.055 0.095 0)
(0.06 0.095 0)
(0.065 0.095 0)
(0.07 0.095 0)
(0.075 0.095 0)
(0.08 0.095 0)
(0.085 0.095 0)
(0.09 0.095 0)
(0.095 0.095 0)
(0.1 0.095 0)
(0 0.1 0)
(0.005 0.1 0)
(0.01 0.1 0)
(0.015 0.1 0)
(0.02 0.1 0)
(0.025 0.1 0)
(0.03 0.1 0)
(0.035 0.1 0)
(0.04 0.1 0)
(0.045 0.1 0)
(0.05 0.1 0)
(0.055 0.1 0)
(0.06 0.1 0)
(0.065 0.1 0)
(0.07 0.1 0)
(0.075 0.1 0)
(0.08 0.1 0)
(0.085 0.1 0)
(0.09 0.1 0)
(0.095 0.1 0)
(0.1 0.1 0)
(0 0 0.01)
(0.005 0 0.01)
(0.01 0 0.01)
(0.015 0 0.01)
(0.02 0 0.01)
(0.025 0 0.01)
(0.03 0 0.01)
(0.035 0 0.01)
(0.04 0 0.01)
(0.045 0 0.01)
(0.05 0 0.01)
(0.055 0 0.01)
(0.06 0 0.01)
(0.065 0 0.01)
(0.07 0 0.01)
(0.075 0 0.01)
(0.08 0 0.01)
(0.085 0 0.01)
(0.09 0 0.01)
(0.095 0 0.01)
(0.1 0 0.01)
(0 0.005 0.01)
(0.005 0.005 0.01)
(0.01 0.005 0.01)
(0.015 0.005 0.01)
(0.02 0.005 0.01)
(0.025 0.005 0.01)
(0.03 0.005 0.01)
(0.035 0.005 0.01)
(0.04 0.005 0.01)
(0.045 0.005 0.01)
(0.05 0.005 0.01)
(0.055 0.005 0.01)
(0.06 0.005 0.01)
(0.065 0.005 0.01)
(0.07 0.005 0.01)
(0.075 0.005 0.01)
(0.08 0.005 0.01)
(0.085 0.005 0.01)
(0.09 0.005 0.01)
(0.095 0.005 0.01)
(0.1 0.005 0.01)
(0 0.01 0.01)
(0.005 0.01 0.01)
(0.01 0.01 0.01)
(0.015 0.01 0.01)
(0.02 0.01 0.01)
(0.025 0.01 0.01)
(0.03 0.01 0.01)
(0.035 0.01 0.01)
(0.04 0.01 0.01)
(0.045 0.01 0.01)
(0.05 0.01 0.01)
(0.055 0.01 0.01)
(0.06 0.01 0.01)
(0.065 0.01 0.01)
(0.07 0.01 0.01)
(0.075 0.01 0.01)
(0.08 0.01 0.01)
(0.085 0.01 0.01)
(0.09 0.01 0.01)
(0.095 0.01 0.01)
(0.1 0.01 0.01)
(0 0.015 0.01)
(0.005 0.015 0.01)
(0.01 0.015 0.01)
(0.015 0.015 0.01)
(0.02 0.015 0.01)
(0.025 0.015 0.01)
(0.03 0.015 0.01)
(0.035 0.015 0.01)
(0.04 0.015 0.01)
(0.045 0.015 0.01)
(0.05 0.015 0.01)
(0.055 0.015 0.01)
(0.06 0.015 0.01)
(0.065 0.015 0.01)
(0.07 0.015 0.01)
(0.075 0.015 0.01)
(0.08 0.015 0.01)
(0.085 0.015 0.01)
(0.09 0.015 0.01)
(0.095 0.015 0.01)
(0.1 0.015 0.01)
(0 0.02 0.01)
(0.005 0.02 0.01)
(0.01 0.02 0.01)
(0.015 0.02 0.01)
(0.02 0.02 0.01)
(0.025 0.02 0.01)
(0.03 0.02 0.01)
(0.035 0.02 0.01)
(0.04 0.02 0.01)
(0.045 0.02 0.01)
(0.05 0.02 0.01)
(0.055 0.02 0.01)
(0.06 0.02 0.01)
(0.065 0.02 0.01)
(0.07 0.02 0.01)
(0.075 0.02 0.01)
(0.08 0.02 0.01)
(0.085 0.02 0.01)
(0.09 0.02 0.01)
(0.095 0.02 0.01)
(0.1 0.02 0.01)
(0 0.025 0.01)
(0.005 0.025 0.01)
(0.01 0.025 0.01)
(0.015 0.025 0.01)
(0.02 0.025 0.01)
(0.025 0.025 0.01)
(0.03 0.025 0.01)
(0.035 0.025 0.01)
(0.04 0.025 0.01)
(0.045 0.025 0.01)
(0.05 0.025 0.01)
(0.055 0.025 0.01)
(0.06 0.025 0.01)
(0.065 0.025 0.01)
(0.07 0.025 0.01)
(0.075 0.025 0.01)
(0.08 0.025 0.01)
(0.085 0.025 0.01)
(0.09 0.025 0.01)
(0.095 0.025 0.01)
(0.1 0.025 0.01)
(0 0.03 0.01)
(0.005 0.03 0.01)
(0.01 0.03 0.01)
(0.015 0.03 0.01)
(0.02 0.03 0.01)
(0.025 0.03 0.01)
(0.03 0.03 0.01)
(0.035 0.03 0.01)
(0.04 0.03 0.01)
(0.045 0.03 0.01)
(0.05 0.03 0.01)
(0.055 0.03 0.01)
(0.06 0.03 0.01)
(0.065 0.03 0.01)
(0.07 0.03 0.01)
(0.075 0.03 0.01)
(0.08 0.03 0.01)
(0.085 0.03 0.01)
(0.09 0.03 0.01)
(0.095 0.03 0.01)
(0.1 0.03 0.01)
(0 0.035 0.01)
(0.005 0.035 0.01)
(0.01 0.035 0.01)
(0.015 0.035 0.01)
(0.02 0.035 0.01)
(0.025 0.035 0.01)
(0.03 0.035 0.01)
(0.035 0.035 0.01)
(0.04 0.035 0.01)
(0.045 0.035 0.01)
(0.05 0.035 0.01)
(0.055 0.035 0.01)
(0.06 0.035 0.01)
(0.065 0.035 0.01)
(0.07 0.035 0.01)
(0.075 0.035 0.01)
(0.08 0.035 0.01)
(0.085 0.035 0.01)
(0.09 0.035 0.01)
(0.095 0.035 0.01)
(0.1 0.035 0.01)
(0 0.04 0.01)
(0.005 0.04 0.01)
(0.01 0.04 0.01)
(0.015 0.04 0.01)
(0.02 0.04 0.01)
(0.025 0.04 0.01)
(0.03 0.04 0.01)
(0.035 0.04 0.01)
(0.04 0.04 0.01)
(0.045 0.04 0.01)
(0.05 0.04 0.01)
(0.055 0.04 0.01)
(0.06 0.04 0.01)
(0.065 0.04 0.01)
(0.07 0.04 0.01)
(0.075 0.04 0.01)
(0.08 0.04 0.01)
(0.085 0.04 0.01)
(0.09 0.04 0.01)
(0.095 0.04 0.01)
(0.1 0.04 0.01)
(0 0.045 0.01)
(0.005 0.045 0.01)
(0.01 0.045 0.01)
(0.015 0.045 0.01)
(0.02 0.045 0.01)
(0.025 0.045 0.01)
(0.03 0.045 0.01)
(0.035 0.045 0.01)
(0.04 0.045 0.01)
(0.045 0.045 0.01)
(0.05 0.045 0.01)
(0.055 0.045 0.01)
(0.06 0.045 0.01)
(0.065 0.045 0.01)
(0.07 0.045 0.01)
(0.075 0.045 0.01)
(0.08 0.045 0.01)
(0.085 0.045 0.01)
(0.09 0.045 0.01)
(0.095 0.045 0.01)
(0.1 0.045 0.01)
(0 0.05 0.01)
(0.005 0.05 0.01)
(0.01 0.05 0.01)
(0.015 0.05 0.01)
(0.02 0.05 0.01)
(0.025 0.05 0.01)
(0.03 0.05 0.01)
(0.035 0.05 0.01)
(0.04 0.05 0.01)
(0.045 0.05 0.01)
(0.05 0.05 0.01)
(0.055 0.05 0.01)
(0.06 0.05 0.01)
(0.065 0.05 0.01)
(0.07 0.05 0.01)
(0.075 0.05 0.01)
(0.08 0.05 0.01)
(0.085 0.05 0.01)
(0.09 0.05 0.01)
(0.095 0.05 0.01)
(0.1 0.05 0.01)
(0 0.055 0.01)
(0.005 0.055 0.01)
(0.01 0.055 0.01)
(0.015 0.055 0.01)
(0.02 0.055 0.01)
(0.025 0.055 0.01)
(0.03 0.055 0.01)
(0.035 0.055 0.01)
(0.04 0.055 0.01)
(0.045 0.055 0.01)
(0.05 0.055 0.01)
(0.055 0.055 0.01)
(0.06 0.055 0.01)
(0.065 0.055 0.01)
(0.07 0.055 0.01)
(0.075 0.055 0.01)
(0.08 0.055 0.01)
(0.085 0.055 0.01)
(0.09 0.055 0.01)
(0.095 0.055 0.01)
(0.1 0.055 0.01)
(0 0.06 0.01)
(0.005 0.06 0.01)
(0.01 0.06 0.01)
(0.015 0.06 0.01)
(0.02 0.06 0.01)
(0.025 0.06 0.01)
(0.03 0.06 0.01)
(0.035 0.06 0.01)
(0.04 0.06 0.01)
(0.045 0.06 0.01)
(0.05 0.06 0.01)
(0.055 0.06 0.01)
(0.06 0.06 0.01)
(0.065 0.06 0.01)
(0.07 0.06 0.01)
(0.075 0.06 0.01)
(0.08 0.06 0.01)
(0.085 0.06 0.01)
(0.09 0.06 0.01)
(0.095 0.06 0.01)
(0.1 0.06 0.01)
(0 0.065 0.01)
(0.005 0.065 0.01)
(0.01 0.065 0.01)
(0.015 0.065 0.01)
(0.02 0.065 0.01)
(0.025 0.065 0.01)
(0.03 0.065 0.01)
(0.035 0.065 0.01)
(0.04 0.065 0.01)
(0.045 0.065 0.01)
(0.05 0.065 0.01)
(0.055 0.065 0.01)
(0.06 0.065 0.01)
(0.065 0.065 0.01)
(0.07 0.065 0.01)
(0.075 0.065 0.01)
(0.08 0.065 0.01)
(0.085 0.065 0.01)
(0.09 0.065 0.01)
(0.095 0.065 0.01)
(0.1 0.065 0.01)
(0 0.07 0.01)
(0.005 0.07 0.01)
(0.01 0.07 0.01)
(0.015 0.07 0.01)
(0.02 0.07 0.01)
(0.025 0.07 0.01)
(0.03 0.07 0.01)
(0.035 0.07 0.01)
(0.04 0.07 0.01)
(0.045 0.07 0.01)
(0.05 0.07 0.01)
(0.055 0.07 0.01)
(0.06 0.07 0.01)
(0.065 0.07 0.01)
(0.07 0.07 0.01)
(0.075 0.07 0.01)
(0.08 0.07 0.01)
(0.085 0.07 0.01)
(0.09 0.07 0.01)
(0.095 0.07 0.01)
(0.1 0.07 0.01)
(0 0.075 0.01)
(0.005 0.075 0.01)
(0.01 0.075 0.01)
(0.015 0.075 0.01)
(0.02 0.075 0.01)
(0.025 0.075 0.01)
(0.03 0.075 0.01)
(0.035 0.075 0.01)
(0.04 0.075 0.01)
(0.045 0.075 0.01)
(0.05 0.075 0.01)
(0.055 0.075 0.01)
(0.06 0.075 0.01)
(0.065 0.075 0.01)
(0.07 0.075 0.01)
(0.075 0.075 0.01)
(0.08 0.075 0.01)
(0.085 0.075 0.01)
(0.09 0.075 0.01)
(0.095 0.075 0.01)
(0.1 0.075 0.01)
(0 0.08 0.01)
(0.005 0.08 0.01)
(0.01 0.08 0.01)
(0.015 0.08 0.01)
(0.02 0.08 0.01)
(0.025 0.08 0.01)
(0.03 0.08 0.01)
(0.035 0.08 0.01)
(0.04 0.08 0.01)
(0.045 0.08 0.01)
(0.05 0.08 0.01)
(0.055 0.08 0.01)
(0.06 0.08 0.01)
(0.065 0.08 0.01)
(0.07 0.08 0.01)
(0.075 0.08 0.01)
(0.08 0.08 0.01)
(0.085 0.08 0.01)
(0.09 0.08 0.01)
(0.095 0.08 0.01)
(0.1 0.08 0.01)
(0 0.085 0.01)
(0.005 0.085 0.01)
(0.01 0.085 0.01)
(0.015 0.085 0.01)
(0.02 0.085 0.01)
(0.025 0.085 0.01)
(0.03 0.085 0.01)
(0.035 0.085 0.01)
(0.04 0.085 0.01)
(0.045 0.085 0.01)
(0.05 0.085 0.01)
(0.055 0.085 0.01)
(0.06 0.085 0.01)
(0.065 0.085 0.01)
(0.07 0.085 0.01)
(0.075 0.085 0.01)
(0.08 0.085 0.01)
(0.085 0.085 0.01)
(0.09 0.085 0.01)
(0.095 0.085 0.01)
(0.1 0.085 0.01)
(0 0.09 0.01)
(0.005 0.09 0.01)
(0.01 0.09 0.01)
(0.015 0.09 0.01)
(0.02 0.09 0.01)
(0.025 0.09 0.01)
(0.03 0.09 0.01)
(0.035 0.09 0.01)
(0.04 0.09 0.01)
(0.045 0.09 0.01)
(0.05 0.09 0.01)
(0.055 0.09 0.01)
(0.06 0.09 0.01)
(0.065 0.09 0.01)
(0.07 0.09 0.01)
(0.075 0.09 0.01)
(0.08 0.09 0.01)
(0.085 0.09 0.01)
(0.09 0.09 0.01)
(0.095 0.09 0.01)
(0.1 0.09 0.01)
(0 0.095 0.01)
(0.005 0.095 0.01)
(0.01 0.095 0.01)
(0.015 0.095 0.01)
(0.02 0.095 0.01)
(0.025 0.095 0.01)
(0.03 0.095 0.01)
(0.035 0.095 0.01)
(0.04 0.095 0.01)
(0.045 0.095 0.01)
(0.05 0.095 0.01)
(0.055 0.095 0.01)
(0.06 0.095 0.01)
(0.065 0.095 0.01)
(0.07 0.095 0.01)
(0.075 0.095 0.01)
(0.08 0.095 0.01)
(0.085 0.095 0.01)
(0.09 0.095 0.01)
(0.095 0.095 0.01)
(0.1 0.095 0.01)
(0 0.1 0.01)
(0.005 0.1 0.01)
(0.01 0.1 0.01)
(0.015 0.1 0.01)
(0.02 0.1 0.01)
(0.025 0.1 0.01)
(0.03 0.1 0.01)
(0.035 0.1 0.01)
(0.04 0.1 0.01)
(0.045 0.1 0.01)
(0.05 0.1 0.01)
(0.055 0.1 0.01)
(0.06 0.1 0.01)
(0.065 0.1 0.01)
(0.07 0.1 0.01)
(0.075 0.1 0.01)
(0.08 0.1 0.01)
(0.085 0.1 0.01)
(0.09 0.1 0.01)
(0.095 0.1 0.01)
(0.1 0.1 0.01)
)
// ************************************************************************* //