diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index d964ac12cc0607feae48d7b56a5c351dac92614b..ef521a257d2238f164acab87bf6e4307d8c8218a 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -803,7 +803,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
if (allGeometry)
{
cellSet cells(mesh, "underdeterminedCells", mesh.nCells()/100);
- if (mesh.checkCellDeterminant(true, &cells, mesh.geometricD()))
+ if (mesh.checkCellDeterminant(true, &cells))
{
noFailedChecks++;
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index 1cd7f37fe88d62c3695395e9944b02fd89e9bb56..55208034a17659f7e023b7ad980f38cb131332ec 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -425,6 +425,10 @@ $(polyMesh)/polyMeshInitMesh.C
$(polyMesh)/polyMeshClear.C
$(polyMesh)/polyMeshUpdate.C
+polyMeshCheck = $(polyMesh)/polyMeshCheck
+$(polyMeshCheck)/polyMeshCheck.C
+$(polyMeshCheck)/polyMeshTools.C
+
primitiveMesh = meshes/primitiveMesh
$(primitiveMesh)/primitiveMesh.C
$(primitiveMesh)/primitiveMeshCellCells.C
@@ -445,9 +449,9 @@ $(primitiveMesh)/primitiveMeshCalcCellShapes.C
primitiveMeshCheck = $(primitiveMesh)/primitiveMeshCheck
$(primitiveMeshCheck)/primitiveMeshCheck.C
-$(primitiveMeshCheck)/primitiveMeshCheckMotion.C
$(primitiveMeshCheck)/primitiveMeshCheckPointNearness.C
$(primitiveMeshCheck)/primitiveMeshCheckEdgeLength.C
+$(primitiveMeshCheck)/primitiveMeshTools.C
primitivePatch = $(primitiveMesh)/primitivePatch
$(primitivePatch)/patchZones.C
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 6388190bf6e5fd10cd9fd38eae5d004d3864bcbd..abe537320f9d0d5b92b2c2589c6047eac269d199 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -1129,7 +1129,7 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
if (debug)
{
// Check mesh motion
- if (primitiveMesh::checkMeshMotion(points_, true))
+ if (checkMeshMotion(points_, true))
{
Info<< "tmp<scalarField> polyMesh::movePoints"
<< "(const pointField&) : "
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 09c046feb840f098c5a4fb80fd3533b13eb84987..e4c189767192ae87ddddc131ea0ea3e0d1068353 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -34,6 +34,7 @@ SourceFiles
polyMeshFromShapeMesh.C
polyMeshIO.C
polyMeshUpdate.C
+ polyMeshCheck.C
\*---------------------------------------------------------------------------*/
@@ -224,7 +225,45 @@ private:
cellList& cells
);
+ // Geometry checks
+ //- Check non-orthogonality
+ bool checkFaceOrthogonality
+ (
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
+ const bool report,
+ const bool detailedReport,
+ labelHashSet* setPtr
+ ) const;
+
+ //- Check face skewness
+ bool checkFaceSkewness
+ (
+ const pointField& points,
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
+ const bool report,
+ const bool detailedReport,
+ labelHashSet* setPtr
+ ) const;
+
+ bool checkEdgeAlignment
+ (
+ const pointField& p,
+ const bool report,
+ const Vector<label>& directions,
+ labelHashSet* setPtr
+ ) const;
+
+ bool checkCellDeterminant
+ (
+ const vectorField& faceAreas,
+ const bool report,
+ labelHashSet* setPtr,
+ const Vector<label>& meshD
+ ) const;
public:
@@ -535,6 +574,48 @@ public:
void removeFiles() const;
+ // Geometric checks. Selectively override primitiveMesh functionality.
+
+ //- Check boundary for closedness
+ virtual bool checkClosedBoundary(const bool report = false) const;
+
+ //- Check non-orthogonality
+ virtual bool checkFaceOrthogonality
+ (
+ const bool report = false,
+ labelHashSet* setPtr = NULL
+ ) const;
+
+ //- Check face skewness
+ virtual bool checkFaceSkewness
+ (
+ const bool report = false,
+ labelHashSet* setPtr = NULL
+ ) const;
+
+ //- Check edge alignment for 1D/2D cases
+ virtual bool checkEdgeAlignment
+ (
+ const bool report,
+ const Vector<label>& directions,
+ labelHashSet* setPtr
+ ) const;
+
+ virtual bool checkCellDeterminant
+ (
+ const bool report,
+ labelHashSet* setPtr
+ ) const;
+
+ //- Check mesh motion for correctness given motion points
+ virtual bool checkMeshMotion
+ (
+ const pointField& newPoints,
+ const bool report = false,
+ const bool detailedReport = false
+ ) const;
+
+
// Helper functions
//- Find the cell, tetFaceI and tetPtI for the given position
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
new file mode 100644
index 0000000000000000000000000000000000000000..748c7da21a63f0b9d25acb21f34d5df839415ee0
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
@@ -0,0 +1,679 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 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 "polyMesh.H"
+#include "polyMeshTools.H"
+#include "unitConversion.H"
+#include "syncTools.H"
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+bool Foam::polyMesh::checkFaceOrthogonality
+(
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
+ const bool report,
+ const bool detailedReport,
+ labelHashSet* setPtr
+) const
+{
+ if (debug)
+ {
+ Info<< "bool polyMesh::checkFaceOrthogonality("
+ << "const bool, labelHashSet*) const: "
+ << "checking mesh non-orthogonality" << endl;
+ }
+
+ const labelList& own = faceOwner();
+ const labelList& nei = faceNeighbour();
+
+ // Calculate orthogonality for all internal and coupled boundary faces
+ // (1 for uncoupled boundary faces)
+ tmp<scalarField> tortho = polyMeshTools::faceOrthogonality
+ (
+ *this,
+ fAreas,
+ cellCtrs
+ );
+ const scalarField& ortho = tortho();
+
+ // Severe nonorthogonality threshold
+ const scalar severeNonorthogonalityThreshold =
+ ::cos(degToRad(primitiveMesh::nonOrthThreshold_));
+
+
+ scalar minDDotS = GREAT;
+ scalar sumDDotS = 0.0;
+ label nSummed = 0;
+ label severeNonOrth = 0;
+ label errorNonOrth = 0;
+
+
+ // Statistics only for internal and masters of coupled faces
+ PackedBoolList isMasterFace(syncTools::getInternalOrMasterFaces(*this));
+
+ forAll(ortho, faceI)
+ {
+ if (ortho[faceI] < severeNonorthogonalityThreshold)
+ {
+ if (ortho[faceI] > SMALL)
+ {
+ if (setPtr)
+ {
+ setPtr->insert(faceI);
+ }
+
+ severeNonOrth++;
+ }
+ else
+ {
+ // Error : non-ortho too large
+ if (setPtr)
+ {
+ setPtr->insert(faceI);
+ }
+ if (detailedReport && errorNonOrth == 0)
+ {
+ // Non-orthogonality greater than 90 deg
+ WarningIn
+ (
+ "polyMesh::checkFaceOrthogonality"
+ "(const pointField&, const bool) const"
+ ) << "Severe non-orthogonality for face "
+ << faceI
+ << " between cells " << own[faceI]
+ << " and " << nei[faceI]
+ << ": Angle = " << radToDeg(::acos(ortho[faceI]))
+ << " deg." << endl;
+ }
+
+ errorNonOrth++;
+ }
+ }
+
+ if (isMasterFace[faceI])
+ {
+ minDDotS = min(minDDotS, ortho[faceI]);
+ sumDDotS += ortho[faceI];
+ nSummed++;
+ }
+ }
+
+ reduce(minDDotS, minOp<scalar>());
+ reduce(sumDDotS, sumOp<scalar>());
+ reduce(nSummed, sumOp<label>());
+ reduce(severeNonOrth, sumOp<label>());
+ reduce(errorNonOrth, sumOp<label>());
+
+ if (debug || report)
+ {
+ if (nSummed > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " Mesh non-orthogonality Max: "
+ << radToDeg(::acos(minDDotS))
+ << " average: " << radToDeg(::acos(sumDDotS/nSummed))
+ << endl;
+ }
+ }
+
+ if (severeNonOrth > 0)
+ {
+ Info<< " *Number of severely non-orthogonal faces: "
+ << severeNonOrth << "." << endl;
+ }
+ }
+
+ if (errorNonOrth > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " ***Number of non-orthogonality errors: "
+ << errorNonOrth << "." << endl;
+ }
+
+ return true;
+ }
+ else
+ {
+ if (debug || report)
+ {
+ Info<< " Non-orthogonality check OK." << endl;
+ }
+
+ return false;
+ }
+}
+
+
+bool Foam::polyMesh::checkFaceSkewness
+(
+ const pointField& points,
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
+ const bool report,
+ const bool detailedReport,
+ labelHashSet* setPtr
+) const
+{
+ if (debug)
+ {
+ Info<< "bool polyMesh::checkFaceSkewnesss("
+ << "const bool, labelHashSet*) const: "
+ << "checking face skewness" << endl;
+ }
+
+ const labelList& own = faceOwner();
+ const labelList& nei = faceNeighbour();
+
+ // Warn if the skew correction vector is more than skewWarning times
+ // larger than the face area vector
+
+ tmp<scalarField> tskew = polyMeshTools::faceSkewness
+ (
+ *this,
+ points,
+ fCtrs,
+ fAreas,
+ cellCtrs
+ );
+ const scalarField& skew = tskew();
+
+ scalar maxSkew = max(skew);
+ label nWarnSkew = 0;
+
+ // Statistics only for all faces except slave coupled faces
+ PackedBoolList isMasterFace(syncTools::getMasterFaces(*this));
+
+ forAll(skew, faceI)
+ {
+ // Check if the skewness vector is greater than the PN vector.
+ // This does not cause trouble but is a good indication of a poor mesh.
+ if (skew[faceI] > skewThreshold_)
+ {
+ if (setPtr)
+ {
+ setPtr->insert(faceI);
+ }
+ if (detailedReport && nWarnSkew == 0)
+ {
+ // Non-orthogonality greater than 90 deg
+ if (isInternalFace(faceI))
+ {
+ WarningIn
+ (
+ "polyMesh::checkFaceSkewnesss"
+ "(const pointField&, const bool) const"
+ ) << "Severe skewness " << skew[faceI]
+ << " for face " << faceI
+ << " between cells " << own[faceI]
+ << " and " << nei[faceI];
+ }
+ else
+ {
+ WarningIn
+ (
+ "polyMesh::checkFaceSkewnesss"
+ "(const pointField&, const bool) const"
+ ) << "Severe skewness " << skew[faceI]
+ << " for boundary face " << faceI
+ << " on cell " << own[faceI];
+ }
+ }
+
+ if (isMasterFace[faceI])
+ {
+ nWarnSkew++;
+ }
+ }
+ }
+
+ reduce(maxSkew, maxOp<scalar>());
+ reduce(nWarnSkew, sumOp<label>());
+
+ if (nWarnSkew > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " ***Max skewness = " << maxSkew
+ << ", " << nWarnSkew << " highly skew faces detected"
+ " which may impair the quality of the results"
+ << endl;
+ }
+
+ return true;
+ }
+ else
+ {
+ if (debug || report)
+ {
+ Info<< " Max skewness = " << maxSkew << " OK." << endl;
+ }
+
+ return false;
+ }
+}
+
+
+// Check 1D/2Dness of edges. Gets passed the non-empty directions and
+// checks all edges in the mesh whether they:
+// - have no component in a non-empty direction or
+// - are only in a singe non-empty direction.
+// Empty direction info is passed in as a vector of labels (synchronised)
+// which are 1 if the direction is non-empty, 0 if it is.
+bool Foam::polyMesh::checkEdgeAlignment
+(
+ const pointField& p,
+ const bool report,
+ const Vector<label>& directions,
+ labelHashSet* setPtr
+) const
+{
+ if (debug)
+ {
+ Info<< "bool polyMesh::checkEdgeAlignment("
+ << "const bool, const Vector<label>&, labelHashSet*) const: "
+ << "checking edge alignment" << endl;
+ }
+
+ label nDirs = 0;
+ for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+ {
+ if (directions[cmpt] == 1)
+ {
+ nDirs++;
+ }
+ else if (directions[cmpt] != 0)
+ {
+ FatalErrorIn
+ (
+ "polyMesh::checkEdgeAlignment"
+ "(const bool, const Vector<label>&, labelHashSet*)"
+ ) << "directions should contain 0 or 1 but is now " << directions
+ << exit(FatalError);
+ }
+ }
+
+ if (nDirs == vector::nComponents)
+ {
+ return false;
+ }
+
+
+ const faceList& fcs = faces();
+
+ EdgeMap<label> edgesInError;
+
+ forAll(fcs, faceI)
+ {
+ const face& f = fcs[faceI];
+
+ forAll(f, fp)
+ {
+ label p0 = f[fp];
+ label p1 = f.nextLabel(fp);
+ if (p0 < p1)
+ {
+ vector d(p[p1]-p[p0]);
+ scalar magD = mag(d);
+
+ if (magD > ROOTVSMALL)
+ {
+ d /= magD;
+
+ // Check how many empty directions are used by the edge.
+ label nEmptyDirs = 0;
+ label nNonEmptyDirs = 0;
+ for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+ {
+ if (mag(d[cmpt]) > 1e-6)
+ {
+ if (directions[cmpt] == 0)
+ {
+ nEmptyDirs++;
+ }
+ else
+ {
+ nNonEmptyDirs++;
+ }
+ }
+ }
+
+ if (nEmptyDirs == 0)
+ {
+ // Purely in ok directions.
+ }
+ else if (nEmptyDirs == 1)
+ {
+ // Ok if purely in empty directions.
+ if (nNonEmptyDirs > 0)
+ {
+ edgesInError.insert(edge(p0, p1), faceI);
+ }
+ }
+ else if (nEmptyDirs > 1)
+ {
+ // Always an error
+ edgesInError.insert(edge(p0, p1), faceI);
+ }
+ }
+ }
+ }
+ }
+
+ label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
+
+ if (nErrorEdges > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " ***Number of edges not aligned with or perpendicular to "
+ << "non-empty directions: " << nErrorEdges << endl;
+ }
+
+ if (setPtr)
+ {
+ setPtr->resize(2*edgesInError.size());
+ forAllConstIter(EdgeMap<label>, edgesInError, iter)
+ {
+ setPtr->insert(iter.key()[0]);
+ setPtr->insert(iter.key()[1]);
+ }
+ }
+
+ return true;
+ }
+ else
+ {
+ if (debug || report)
+ {
+ Info<< " All edges aligned with or perpendicular to "
+ << "non-empty directions." << endl;
+ }
+ return false;
+ }
+}
+
+
+bool Foam::polyMesh::checkCellDeterminant
+(
+ const vectorField& faceAreas,
+ const bool report,
+ labelHashSet* setPtr,
+ const Vector<label>& meshD
+) const
+{
+ if (debug)
+ {
+ Info<< "bool polyMesh::checkCellDeterminant(const bool"
+ << ", labelHashSet*) const: "
+ << "checking for under-determined cells" << endl;
+ }
+
+ // 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 = primitiveMeshTools::cellDeterminant
+ (
+ *this,
+ meshD,
+ faceAreas,
+ syncTools::getInternalOrCoupledFaces(*this)
+ );
+ scalarField& cellDeterminant = tcellDeterminant();
+
+
+ label nErrorCells = 0;
+ scalar minDet = min(cellDeterminant);
+ scalar sumDet = sum(cellDeterminant);
+
+ forAll (cellDeterminant, cellI)
+ {
+ if (cellDeterminant[cellI] < 1e-3)
+ {
+ if (setPtr)
+ {
+ setPtr->insert(cellI);
+ }
+
+ nErrorCells++;
+ }
+ }
+
+ reduce(nErrorCells, sumOp<label>());
+ reduce(minDet, minOp<scalar>());
+ reduce(sumDet, sumOp<scalar>());
+ label nSummed = returnReduce(cellDeterminant.size(), sumOp<label>());
+
+ if (debug || report)
+ {
+ if (nSummed > 0)
+ {
+ Info<< " Cell determinant (wellposedness) : minimum: " << minDet
+ << " average: " << sumDet/nSummed
+ << endl;
+ }
+ }
+
+ if (nErrorCells > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " ***Cells with small determinant found, number of cells: "
+ << nErrorCells << endl;
+ }
+
+ return true;
+ }
+ else
+ {
+ if (debug || report)
+ {
+ Info<< " Cell determinant check OK." << endl;
+ }
+
+ return false;
+ }
+
+ return false;
+}
+
+
+bool Foam::polyMesh::checkClosedBoundary(const bool report) const
+{
+ return primitiveMesh::checkClosedBoundary
+ (
+ faceAreas(),
+ report,
+ syncTools::getInternalOrCoupledFaces(*this)
+ );
+}
+
+
+bool Foam::polyMesh::checkFaceOrthogonality
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ return checkFaceOrthogonality
+ (
+ faceAreas(),
+ cellCentres(),
+ report,
+ false, // detailedReport
+ setPtr
+ );
+}
+
+
+bool Foam::polyMesh::checkFaceSkewness
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ return checkFaceSkewness
+ (
+ points(),
+ faceCentres(),
+ faceAreas(),
+ cellCentres(),
+ report,
+ false, // detailedReport
+ setPtr
+ );
+}
+
+
+bool Foam::polyMesh::checkEdgeAlignment
+(
+ const bool report,
+ const Vector<label>& directions,
+ labelHashSet* setPtr
+) const
+{
+ return checkEdgeAlignment
+ (
+ points(),
+ report,
+ directions,
+ setPtr
+ );
+}
+
+
+bool Foam::polyMesh::checkCellDeterminant
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ return checkCellDeterminant
+ (
+ faceAreas(),
+ report,
+ setPtr,
+ geometricD()
+ );
+}
+
+
+bool Foam::polyMesh::checkMeshMotion
+(
+ const pointField& newPoints,
+ const bool report,
+ const bool detailedReport
+) const
+{
+ if (debug || report)
+ {
+ Pout<< "bool polyMesh::checkMeshMotion("
+ << "const pointField&, const bool, const bool) const: "
+ << "checking mesh motion" << endl;
+ }
+
+ vectorField fCtrs(nFaces());
+ vectorField fAreas(nFaces());
+
+ makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
+
+ // Check cell volumes and calculate new cell centres
+ vectorField cellCtrs(nCells());
+ scalarField cellVols(nCells());
+
+ makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
+
+ // Check cell volumes
+ bool error = checkCellVolumes
+ (
+ cellVols, // vols
+ report, // report
+ detailedReport, // detailedReport
+ NULL // setPtr
+ );
+
+
+ // Check face areas
+ bool areaError = checkFaceAreas
+ (
+ faceAreas(),
+ report, // report
+ detailedReport, // detailedReport,
+ NULL // setPtr
+ );
+ error = error || areaError;
+
+
+ // Check pyramid volumes
+ bool pyrVolError = checkFacePyramids
+ (
+ newPoints,
+ cellCtrs,
+ report, // report,
+ detailedReport, // detailedReport,
+ -SMALL, // minPyrVol
+ NULL // setPtr
+ );
+ error = error || pyrVolError;
+
+
+ // Check face non-orthogonality
+ bool nonOrthoError = checkFaceOrthogonality
+ (
+ fAreas,
+ cellCtrs,
+ report, // report
+ detailedReport, // detailedReport
+ NULL // setPtr
+ );
+ error = error || nonOrthoError;
+
+
+ if (!error && (debug || report))
+ {
+ Pout<< "Mesh motion check OK." << endl;
+ }
+
+ return error;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
new file mode 100644
index 0000000000000000000000000000000000000000..fc0739e79faf42d048b133d79b2a97f572b99388
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
@@ -0,0 +1,190 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 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::polyMeshTools::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();
+
+ // Internal faces
+ forAll(nei, faceI)
+ {
+ ortho[faceI] = primitiveMeshTools::faceOrthogonality
+ (
+ cc[own[faceI]],
+ cc[nei[faceI]],
+ areas[faceI]
+ );
+ }
+
+
+ // Coupled faces
+
+ pointField neighbourCc;
+ syncTools::swapBoundaryCellList(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] = primitiveMeshTools::faceOrthogonality
+ (
+ cc[own[faceI]],
+ neighbourCc[bFaceI],
+ areas[faceI]
+ );
+ }
+ }
+ }
+
+ return tortho;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::polyMeshTools::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 faceList& fcs = mesh.faces();
+ const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+ tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
+ scalarField& skew = tskew();
+
+ forAll(nei, faceI)
+ {
+ skew[faceI] = primitiveMeshTools::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::swapBoundaryCellList(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] = primitiveMeshTools::faceSkewness
+ (
+ mesh,
+ p,
+ fCtrs,
+ fAreas,
+
+ faceI,
+ cellCtrs[own[faceI]],
+ neighbourCc[bFaceI]
+ );
+ }
+ }
+ else
+ {
+ forAll(pp, i)
+ {
+ label faceI = pp.start() + i;
+
+ vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
+
+ vector normal = fAreas[faceI];
+ normal /= mag(normal) + VSMALL;
+ vector d = normal*(normal & Cpf);
+
+
+ // Skewness vector
+ vector sv =
+ Cpf
+ - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
+ vector svHat = sv/(mag(sv) + VSMALL);
+
+ // 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) + VSMALL;
+ const face& f = fcs[faceI];
+ forAll(f, pi)
+ {
+ fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
+ }
+
+ // Normalised skewness
+ skew[faceI] = mag(sv)/fd;
+ }
+ }
+ }
+
+ return tskew;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
new file mode 100644
index 0000000000000000000000000000000000000000..53c758c3be7967f5a914b39e7e0c12a50bbaf86c
--- /dev/null
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 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::polyMeshTools
+
+Description
+ Collection of static 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 polyMeshTools Declaration
+\*---------------------------------------------------------------------------*/
+
+class polyMeshTools
+:
+ public primitiveMeshTools
+{
+
+public:
+
+ //- Generate orthogonality field. (1 for fully orthogonal, < 1 for
+ // non-orthogonal)
+ static tmp<scalarField> faceOrthogonality
+ (
+ const polyMesh& mesh,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs
+ );
+
+// static tmp<scalarField> faceOrthogonality(const polyMesh& mesh)
+// {
+// return faceOrthogonality
+// (
+// mesh,
+// mesh.faceAreas(),
+// mesh.cellCentres()
+// );
+// }
+
+ //- Generate skewness field
+ static tmp<scalarField> faceSkewness
+ (
+ const polyMesh& mesh,
+ const pointField& points,
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs
+ );
+
+// static tmp<scalarField> faceSkewness(const polyMesh& mesh)
+// {
+// return faceSkewness
+// (
+// mesh,
+// mesh.points(),
+// mesh.faceCentres(),
+// mesh.faceAreas(),
+// mesh.cellCentres()
+// );
+// }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index 8c647d8454c565f2432ebb877f9c73a8c703c6ae..c4ab139ca3f99f4b8a9366e9ede3bc7444c76853 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -43,9 +43,6 @@ SourceFiles
primitiveMeshEdges.C
primitiveMeshCellCentresAndVols.C
primitiveMeshFaceCentresAndAreas.C
- primitiveMeshEdgeVectors.C
- primitiveMeshCheck.C
- primitiveMeshCheckMotion.C
primitiveMeshFindCell.C
\*---------------------------------------------------------------------------*/
@@ -69,6 +66,8 @@ SourceFiles
namespace Foam
{
+class PackedBoolList;
+
/*---------------------------------------------------------------------------*\
Class primitiveMesh Declaration
\*---------------------------------------------------------------------------*/
@@ -225,6 +224,28 @@ class primitiveMesh
const labelList&
);
+protected:
+
+ // Static data members
+
+ //- Static data to control mesh checking
+
+ //- Cell closedness warning threshold
+ // set as the fraction of un-closed area to closed area
+ static scalar closedThreshold_;
+
+ //- Aspect ratio warning threshold
+ static scalar aspectThreshold_;
+
+ //- Non-orthogonality warning threshold in deg
+ static scalar nonOrthThreshold_;
+
+ //- Skewness warning threshold
+ static scalar skewThreshold_;
+
+ //- Threshold where faces are considered coplanar
+ static scalar planarCosAngle_;
+
// Geometrical calculations
@@ -251,7 +272,7 @@ class primitiveMesh
void calcEdgeVectors() const;
- // Helper functions for mesh checking
+ // Mesh checking
//- Check if all points on face are shared with another face.
bool checkDuplicateFaces
@@ -270,29 +291,104 @@ class primitiveMesh
labelHashSet*
) const;
+ //- Check boundary for closedness
+ bool checkClosedBoundary
+ (
+ const vectorField&,
+ const bool,
+ const PackedBoolList&
+ ) const;
- // Static data members
+ //- Check cells for closedness
+ bool checkClosedCells
+ (
+ const vectorField& faceAreas,
+ const scalarField& cellVolumes,
+ const bool report,
+ labelHashSet* setPtr,
+ labelHashSet* aspectSetPtr,
+ const Vector<label>& meshD
+ ) const;
- //- Static data to control mesh checking
+ //- Check for negative face areas
+ bool checkFaceAreas
+ (
+ const vectorField& faceAreas,
+ const bool report,
+ const bool detailedReport,
+ labelHashSet* setPtr
+ ) const;
- //- Cell closedness warning threshold
- // set as the fraction of un-closed area to closed area
- static scalar closedThreshold_;
+ //- Check for negative cell volumes
+ bool checkCellVolumes
+ (
+ const scalarField& vols,
+ const bool report,
+ const bool detailedReport,
+ labelHashSet* setPtr
+ ) const;
- //- Aspect ratio warning threshold
- static scalar aspectThreshold_;
+ //- Check for non-orthogonality
+ bool checkFaceOrthogonality
+ (
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
+ const bool report,
+ labelHashSet* setPtr
+ ) const;
- //- Non-orthogonality warning threshold in deg
- static scalar nonOrthThreshold_;
+ //- Check face pyramid volume
+ bool checkFacePyramids
+ (
+ const pointField& points,
+ const vectorField& ctrs,
+ const bool report,
+ const bool detailedReport,
+ const scalar minPyrVol,
+ labelHashSet* setPtr
+ ) const;
- //- Skewness warning threshold
- static scalar skewThreshold_;
+ //- Check face skewness
+ bool checkFaceSkewness
+ (
+ const pointField& points,
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
+ const bool report,
+ labelHashSet* setPtr
+ ) const;
- //- Threshold where faces are considered coplanar
- static scalar planarCosAngle_;
+ //- Check face angles
+ bool checkFaceAngles
+ (
+ const pointField& points,
+ const vectorField& faceAreas,
+ const bool report,
+ const scalar maxDeg,
+ labelHashSet* setPtr
+ ) const;
+ //- Check face warpage
+ bool checkFaceFlatness
+ (
+ const pointField& points,
+ const vectorField& faceCentres,
+ const vectorField& faceAreas,
+ const bool report,
+ const scalar warnFlatness,
+ labelHashSet* setPtr
+ ) const;
+
+ //- Check for concave cells by the planes of faces
+ bool checkConcaveCells
+ (
+ const vectorField& fAreas,
+ const pointField& fCentres,
+ const bool report,
+ labelHashSet* setPtr
+ ) const;
-protected:
//- Construct null
primitiveMesh();
@@ -502,29 +598,36 @@ public:
// Topological checks
+ //- Check face ordering
+ virtual bool checkUpperTriangular
+ (
+ const bool report = false,
+ labelHashSet* setPtr = NULL
+ ) const;
+
//- Check cell zip-up
- bool checkCellsZipUp
+ virtual bool checkCellsZipUp
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
//- Check uniqueness of face vertices
- bool checkFaceVertices
+ virtual bool checkFaceVertices
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
- //- Check face-face connectivity
- bool checkFaceFaces
+ //- Check for unused points
+ virtual bool checkPoints
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
- //- Check face ordering
- bool checkUpperTriangular
+ //- Check face-face connectivity
+ virtual bool checkFaceFaces
(
const bool report = false,
labelHashSet* setPtr = NULL
@@ -534,10 +637,11 @@ public:
// Geometric checks
//- Check boundary for closedness
- bool checkClosedBoundary(const bool report = false) const;
+ virtual bool checkClosedBoundary(const bool report = false)
+ const;
//- Check cells for closedness
- bool checkClosedCells
+ virtual bool checkClosedCells
(
const bool report = false,
labelHashSet* setPtr = NULL,
@@ -546,28 +650,28 @@ public:
) const;
//- Check for negative face areas
- bool checkFaceAreas
+ virtual bool checkFaceAreas
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
//- Check for negative cell volumes
- bool checkCellVolumes
+ virtual bool checkCellVolumes
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
//- Check for non-orthogonality
- bool checkFaceOrthogonality
+ virtual bool checkFaceOrthogonality
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
//- Check face pyramid volume
- bool checkFacePyramids
+ virtual bool checkFacePyramids
(
const bool report = false,
const scalar minPyrVol = -SMALL,
@@ -575,14 +679,14 @@ public:
) const;
//- Check face skewness
- bool checkFaceSkewness
+ virtual bool checkFaceSkewness
(
const bool report = false,
labelHashSet* setPtr = NULL
) const;
//- Check face angles
- bool checkFaceAngles
+ virtual bool checkFaceAngles
(
const bool report = false,
const scalar maxSin = 10, // In degrees
@@ -592,31 +696,16 @@ public:
//- Check face warpage: decompose face and check ratio between
// magnitude of sum of triangle areas and sum of magnitude of
// triangle areas.
- bool checkFaceFlatness
+ virtual bool checkFaceFlatness
(
const bool report,
const scalar warnFlatness, // When to include in set.
labelHashSet* setPtr
) const;
- //- Check edge alignment for 1D/2D cases
- bool checkEdgeAlignment
- (
- const bool report,
- const Vector<label>& directions,
- labelHashSet* setPtr = NULL
- ) const;
-
- //- Check for unused points
- bool checkPoints
- (
- const bool report = false,
- labelHashSet* setPtr = NULL
- ) const;
-
//- Check for point-point-nearness,
// e.g. colocated points which may be part of baffles.
- bool checkPointNearness
+ virtual bool checkPointNearness
(
const bool report,
const scalar reportDistSqr,
@@ -624,23 +713,15 @@ public:
) const;
//- Check edge length
- bool checkEdgeLength
+ virtual bool checkEdgeLength
(
const bool report,
const scalar minLenSqr,
labelHashSet* setPtr = NULL
) const;
- //- Check cell determinant
- bool checkCellDeterminant
- (
- const bool report = false,
- labelHashSet* setPtr = NULL,
- const Vector<label>& solutionD = Vector<label>::one
- ) const;
-
//- Check for concave cells by the planes of faces
- bool checkConcaveCells
+ virtual bool checkConcaveCells
(
const bool report = false,
labelHashSet* setPtr = NULL
@@ -649,22 +730,14 @@ public:
//- Check mesh topology for correctness.
// Returns false for no error.
- bool checkTopology(const bool report = false) const;
+ virtual bool checkTopology(const bool report = false) const;
//- Check mesh geometry (& implicitly topology) for correctness.
// Returns false for no error.
- bool checkGeometry(const bool report = false) const;
+ virtual bool checkGeometry(const bool report = false) const;
//- Check mesh for correctness. Returns false for no error.
- bool checkMesh(const bool report = false) const;
-
- //- Check mesh motion for correctness given motion points
- bool checkMeshMotion
- (
- const pointField& newPoints,
- const bool report = false
- ) const;
-
+ virtual bool checkMesh(const bool report = false) const;
//- Set the closedness ratio warning threshold
static scalar setClosedThreshold(const scalar);
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 76fa1aa5d1d38f2c3840c7975fbcc9f03293f206..d805608e4480c11a75c1605a6d5a7b2fe2fb3877 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -29,6 +29,7 @@ License
#include "unitConversion.H"
#include "SortableList.H"
#include "EdgeMap.H"
+#include "primitiveMeshTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -41,7 +42,12 @@ Foam::scalar Foam::primitiveMesh::planarCosAngle_ = 1.0e-6;
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
+bool Foam::primitiveMesh::checkClosedBoundary
+(
+ const vectorField& areas,
+ const bool report,
+ const PackedBoolList& internalOrCoupledFaces
+) const
{
if (debug)
{
@@ -56,12 +62,13 @@ bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
vector sumClosed(vector::zero);
scalar sumMagClosedBoundary = 0;
- const vectorField& areas = faceAreas();
-
for (label faceI = nInternalFaces(); faceI < areas.size(); faceI++)
{
- sumClosed += areas[faceI];
- sumMagClosedBoundary += mag(areas[faceI]);
+ if (internalOrCoupledFaces.size() && !internalOrCoupledFaces[faceI])
+ {
+ sumClosed += areas[faceI];
+ sumMagClosedBoundary += mag(areas[faceI]);
+ }
}
reduce(sumClosed, sumOp<vector>());
@@ -95,6 +102,8 @@ bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
bool Foam::primitiveMesh::checkClosedCells
(
+ const vectorField& faceAreas,
+ const scalarField& cellVolumes,
const bool report,
labelHashSet* setPtr,
labelHashSet* aspectSetPtr,
@@ -140,66 +149,28 @@ bool Foam::primitiveMesh::checkClosedCells
return true;
}
- // Loop through cell faces and sum up the face area vectors for each cell.
- // This should be zero in all vector components
-
- vectorField sumClosed(nCells(), vector::zero);
- vectorField sumMagClosed(nCells(), vector::zero);
- const labelList& own = faceOwner();
- const labelList& nei = faceNeighbour();
- const vectorField& areas = faceAreas();
-
- 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]);
- }
+ scalarField openness;
+ scalarField aspectRatio;
+ primitiveMeshTools::cellClosedness
+ (
+ *this,
+ meshD,
+ faceAreas,
+ cellVolumes,
+ openness,
+ aspectRatio
+ );
label nOpen = 0;
- scalar maxOpennessCell = 0;
-
+ scalar maxOpennessCell = max(openness);
label nAspect = 0;
- scalar maxAspectRatio = 0;
-
- const scalarField& vols = cellVolumes();
-
- label nDims = 0;
- for (direction dir = 0; dir < vector::nComponents; dir++)
- {
- if (meshD[dir] == 1)
- {
- nDims++;
- }
- }
-
+ scalar maxAspectRatio = max(aspectRatio);
// Check the sums
- forAll(sumClosed, cellI)
+ forAll(openness, cellI)
{
- scalar maxOpenness = 0;
-
- for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
- {
- maxOpenness = max
- (
- maxOpenness,
- mag(sumClosed[cellI][cmpt])
- /(sumMagClosed[cellI][cmpt] + VSMALL)
- );
- }
-
- maxOpennessCell = max(maxOpennessCell, maxOpenness);
-
- if (maxOpenness > closedThreshold_)
+ if (openness[cellI] > closedThreshold_)
{
if (setPtr)
{
@@ -209,32 +180,7 @@ bool Foam::primitiveMesh::checkClosedCells
nOpen++;
}
- // 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 + VSMALL);
- if (nDims == 3)
- {
- aspectRatio = max
- (
- aspectRatio,
- 1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
- );
- }
-
- maxAspectRatio = max(maxAspectRatio, aspectRatio);
-
- if (aspectRatio > aspectThreshold_)
+ if (aspectRatio[cellI] > aspectThreshold_)
{
if (aspectSetPtr)
{
@@ -290,7 +236,9 @@ bool Foam::primitiveMesh::checkClosedCells
bool Foam::primitiveMesh::checkFaceAreas
(
+ const vectorField& faceAreas,
const bool report,
+ const bool detailedReport,
labelHashSet* setPtr
) const
{
@@ -301,7 +249,7 @@ bool Foam::primitiveMesh::checkFaceAreas
<< "checking face area magnitudes" << endl;
}
- const scalarField magFaceAreas(mag(faceAreas()));
+ const scalarField magFaceAreas(mag(faceAreas));
scalar minArea = GREAT;
scalar maxArea = -GREAT;
@@ -314,6 +262,25 @@ bool Foam::primitiveMesh::checkFaceAreas
{
setPtr->insert(faceI);
}
+ if (detailedReport)
+ {
+ if (isInternalFace(faceI))
+ {
+ Pout<< "Zero or negative face area detected for "
+ << "internal face "<< faceI << " between cells "
+ << faceOwner()[faceI] << " and "
+ << faceNeighbour()[faceI]
+ << ". Face area magnitude = " << magFaceAreas[faceI]
+ << endl;
+ }
+ else
+ {
+ Pout<< "Zero or negative face area detected for "
+ << "boundary face " << faceI << " next to cell "
+ << faceOwner()[faceI] << ". Face area magnitude = "
+ << magFaceAreas[faceI] << endl;
+ }
+ }
}
minArea = min(minArea, magFaceAreas[faceI]);
@@ -349,7 +316,9 @@ bool Foam::primitiveMesh::checkFaceAreas
bool Foam::primitiveMesh::checkCellVolumes
(
+ const scalarField& vols,
const bool report,
+ const bool detailedReport,
labelHashSet* setPtr
) const
{
@@ -360,8 +329,6 @@ bool Foam::primitiveMesh::checkCellVolumes
<< "checking cell volumes" << endl;
}
- const scalarField& vols = cellVolumes();
-
scalar minVolume = GREAT;
scalar maxVolume = -GREAT;
@@ -375,6 +342,11 @@ bool Foam::primitiveMesh::checkCellVolumes
{
setPtr->insert(cellI);
}
+ if (detailedReport)
+ {
+ Pout<< "Zero or negative cell volume detected for cell "
+ << cellI << ". Volume = " << vols[cellI] << endl;
+ }
nNegVolCells++;
}
@@ -416,6 +388,8 @@ bool Foam::primitiveMesh::checkCellVolumes
bool Foam::primitiveMesh::checkFaceOrthogonality
(
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
const bool report,
labelHashSet* setPtr
) const
@@ -427,35 +401,33 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
<< "checking mesh non-orthogonality" << endl;
}
- // for all internal faces check that the d dot S product is positive
- const vectorField& centres = cellCentres();
- const vectorField& areas = faceAreas();
- const labelList& own = faceOwner();
- const labelList& nei = faceNeighbour();
+ tmp<scalarField> tortho = primitiveMeshTools::faceOrthogonality
+ (
+ *this,
+ fAreas,
+ cellCtrs
+ );
+ const scalarField& ortho = tortho();
// Severe nonorthogonality threshold
const scalar severeNonorthogonalityThreshold =
::cos(degToRad(nonOrthThreshold_));
- scalar minDDotS = GREAT;
+ scalar minDDotS = min(ortho);
- scalar sumDDotS = 0;
+ scalar sumDDotS = sum(ortho);
label severeNonOrth = 0;
label errorNonOrth = 0;
- forAll(nei, faceI)
- {
- vector d = centres[nei[faceI]] - centres[own[faceI]];
- const vector& s = areas[faceI];
-
- scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
- if (dDotS < severeNonorthogonalityThreshold)
+ forAll(ortho, faceI)
+ {
+ if (ortho[faceI] < severeNonorthogonalityThreshold)
{
- if (dDotS > SMALL)
+ if (ortho[faceI] > SMALL)
{
if (setPtr)
{
@@ -474,13 +446,6 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
errorNonOrth++;
}
}
-
- if (dDotS < minDDotS)
- {
- minDDotS = dDotS;
- }
-
- sumDDotS += dDotS;
}
reduce(minDDotS, minOp<scalar>());
@@ -490,7 +455,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
if (debug || report)
{
- label neiSize = nei.size();
+ label neiSize = ortho.size();
reduce(neiSize, sumOp<label>());
if (neiSize > 0)
@@ -535,7 +500,10 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
bool Foam::primitiveMesh::checkFacePyramids
(
+ const pointField& points,
+ const vectorField& ctrs,
const bool report,
+ const bool detailedReport,
const scalar minPyrVol,
labelHashSet* setPtr
) const
@@ -547,46 +515,63 @@ bool Foam::primitiveMesh::checkFacePyramids
<< "checking face orientation" << endl;
}
- // check whether face area vector points to the cell with higher label
- const vectorField& ctrs = cellCentres();
-
const labelList& own = faceOwner();
const labelList& nei = faceNeighbour();
-
const faceList& f = faces();
- const pointField& p = points();
+
+ scalarField ownPyrVol;
+ scalarField neiPyrVol;
+ primitiveMeshTools::facePyramidVolume
+ (
+ *this,
+ points,
+ ctrs,
+ ownPyrVol,
+ neiPyrVol
+ );
+
label nErrorPyrs = 0;
- forAll(f, faceI)
+ forAll(ownPyrVol, faceI)
{
- // Create the owner pyramid - it will have negative volume
- scalar pyrVol = pyramidPointFaceRef(f[faceI], ctrs[own[faceI]]).mag(p);
-
- if (pyrVol > -minPyrVol)
+ if (ownPyrVol[faceI] < minPyrVol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
+ if (detailedReport)
+ {
+ Pout<< "Negative pyramid volume: " << ownPyrVol[faceI]
+ << " for face " << faceI << " " << f[faceI]
+ << " and owner cell: " << own[faceI] << endl
+ << "Owner cell vertex labels: "
+ << cells()[own[faceI]].labels(faces())
+ << endl;
+ }
nErrorPyrs++;
}
if (isInternalFace(faceI))
{
- // Create the neighbour pyramid - it will have positive volume
- scalar pyrVol =
- pyramidPointFaceRef(f[faceI], ctrs[nei[faceI]]).mag(p);
-
- if (pyrVol < minPyrVol)
+ if (neiPyrVol[faceI] < minPyrVol)
{
if (setPtr)
{
setPtr->insert(faceI);
}
-
+ if (detailedReport)
+ {
+ Pout<< "Negative pyramid volume: " << neiPyrVol[faceI]
+ << " for face " << faceI << " " << f[faceI]
+ << " and neighbour cell: " << nei[faceI] << nl
+ << "Neighbour cell vertex labels: "
+ << cells()[nei[faceI]].labels(faces())
+ << endl;
+ }
nErrorPyrs++;
}
}
@@ -619,6 +604,10 @@ bool Foam::primitiveMesh::checkFacePyramids
bool Foam::primitiveMesh::checkFaceSkewness
(
+ const pointField& points,
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
const bool report,
labelHashSet* setPtr
) const
@@ -633,93 +622,24 @@ bool Foam::primitiveMesh::checkFaceSkewness
// Warn if the skew correction vector is more than skewWarning times
// larger than the face area vector
- const pointField& p = points();
- const faceList& fcs = faces();
-
- const labelList& own = faceOwner();
- const labelList& nei = faceNeighbour();
- const vectorField& cellCtrs = cellCentres();
- const vectorField& faceCtrs = faceCentres();
- const vectorField& fAreas = faceAreas();
-
- scalar maxSkew = 0;
+ tmp<scalarField> tskewness = primitiveMeshTools::faceSkewness
+ (
+ *this,
+ points,
+ fCtrs,
+ fAreas,
+ cellCtrs
+ );
+ const scalarField& skewness = tskewness();
+
+ scalar maxSkew = max(skewness);
label nWarnSkew = 0;
- forAll(nei, faceI)
- {
- vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
- vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
-
- // Skewness vector
- vector sv =
- Cpf - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*d;
- vector svHat = sv/(mag(sv) + VSMALL);
-
- // 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) + VSMALL;
- const face& f = fcs[faceI];
- forAll(f, pi)
- {
- fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
- }
-
- // Normalised skewness
- scalar skewness = mag(sv)/fd;
-
- // Check if the skewness vector is greater than the PN vector.
- // This does not cause trouble but is a good indication of a poor mesh.
- if (skewness > skewThreshold_)
- {
- if (setPtr)
- {
- setPtr->insert(faceI);
- }
-
- nWarnSkew++;
- }
-
- if (skewness > maxSkew)
- {
- maxSkew = skewness;
- }
- }
-
-
- // Boundary faces: consider them to have only skewness error.
- // (i.e. treat as if mirror cell on other side)
-
- for (label faceI = nInternalFaces(); faceI < nFaces(); faceI++)
+ forAll(skewness, faceI)
{
- vector Cpf = faceCtrs[faceI] - cellCtrs[own[faceI]];
-
- vector normal = fAreas[faceI];
- normal /= mag(normal) + VSMALL;
- vector d = normal*(normal & Cpf);
-
-
- // Skewness vector
- vector sv =
- Cpf - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
- vector svHat = sv/(mag(sv) + VSMALL);
-
- // 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) + VSMALL;
- const face& f = fcs[faceI];
- forAll(f, pi)
- {
- fd = max(fd, mag(svHat & (p[f[pi]] - faceCtrs[faceI])));
- }
-
- // Normalised skewness
- scalar skewness = mag(sv)/fd;
-
// Check if the skewness vector is greater than the PN vector.
// This does not cause trouble but is a good indication of a poor mesh.
- if (skewness > skewThreshold_)
+ if (skewness[faceI] > skewThreshold_)
{
if (setPtr)
{
@@ -728,14 +648,8 @@ bool Foam::primitiveMesh::checkFaceSkewness
nWarnSkew++;
}
-
- if (skewness > maxSkew)
- {
- maxSkew = skewness;
- }
}
-
reduce(maxSkew, maxOp<scalar>());
reduce(nWarnSkew, sumOp<label>());
@@ -763,86 +677,14 @@ bool Foam::primitiveMesh::checkFaceSkewness
}
-bool Foam::primitiveMesh::checkPoints
-(
- const bool report,
- labelHashSet* setPtr
-) const
-{
- if (debug)
- {
- Info<< "bool primitiveMesh::checkPoints"
- << "(const bool, labelHashSet*) const: "
- << "checking points" << endl;
- }
-
- label nFaceErrors = 0;
- label nCellErrors = 0;
-
- const labelListList& pf = pointFaces();
-
- forAll(pf, pointI)
- {
- if (pf[pointI].empty())
- {
- if (setPtr)
- {
- setPtr->insert(pointI);
- }
-
- nFaceErrors++;
- }
- }
-
-
- forAll(pf, pointI)
- {
- const labelList& pc = pointCells(pointI);
-
- if (pc.empty())
- {
- if (setPtr)
- {
- setPtr->insert(pointI);
- }
-
- nCellErrors++;
- }
- }
-
- reduce(nFaceErrors, sumOp<label>());
- reduce(nCellErrors, sumOp<label>());
-
- if (nFaceErrors > 0 || nCellErrors > 0)
- {
- if (debug || report)
- {
- Info<< " ***Unused points found in the mesh, "
- "number unused by faces: " << nFaceErrors
- << " number unused by cells: " << nCellErrors
- << endl;
- }
-
- return true;
- }
- else
- {
- if (debug || report)
- {
- Info<< " Point usage OK." << endl;
- }
-
- return false;
- }
-}
-
-
// Check convexity of angles in a face. Allow a slight non-convexity.
// E.g. maxDeg = 10 allows for angles < 190 (or 10 degrees concavity)
// (if truly concave and points not visible from face centre the face-pyramid
// check in checkMesh will fail)
bool Foam::primitiveMesh::checkFaceAngles
(
+ const pointField& points,
+ const vectorField& faceAreas,
const bool report,
const scalar maxDeg,
labelHashSet* setPtr
@@ -867,71 +709,30 @@ bool Foam::primitiveMesh::checkFaceAngles
const scalar maxSin = Foam::sin(degToRad(maxDeg));
- const pointField& p = points();
- const faceList& fcs = faces();
- vectorField faceNormals(faceAreas());
- faceNormals /= mag(faceNormals) + VSMALL;
- scalar maxEdgeSin = 0.0;
+ tmp<scalarField> tfaceAngles = primitiveMeshTools::faceConcavity
+ (
+ maxSin,
+ *this,
+ points,
+ faceAreas
+ );
+ const scalarField& faceAngles = tfaceAngles();
- label nConcave = 0;
+ scalar maxEdgeSin = max(faceAngles);
- label errorFaceI = -1;
+ label nConcave = 0;
- forAll(fcs, faceI)
+ forAll(faceAngles, 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 + VSMALL;
-
- forAll(f, fp0)
+ if (faceAngles[faceI] > SMALL)
{
- // Get vertex after fp
- label fp1 = f.fcIndex(fp0);
-
- // Normalized vector between two consecutive points
- vector e10(p[f[fp1]] - p[f[fp0]]);
- scalar magE10 = mag(e10);
- e10 /= magE10 + VSMALL;
+ nConcave++;
- if (magEPrev > SMALL && magE10 > SMALL)
+ if (setPtr)
{
- 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)
- {
- if (faceI != errorFaceI)
- {
- // Count only one error per face.
- errorFaceI = faceI;
- nConcave++;
- }
-
- if (setPtr)
- {
- setPtr->insert(faceI);
- }
-
- maxEdgeSin = max(maxEdgeSin, magEdgeNormal);
- }
- }
+ setPtr->insert(faceI);
}
-
- ePrev = e10;
- magEPrev = magE10;
}
}
@@ -965,11 +766,11 @@ bool Foam::primitiveMesh::checkFaceAngles
}
-// Check warpage of faces. Is calculated as the difference between areas of
-// individual triangles and the overall area of the face (which ifself is
-// is the average of the areas of the individual triangles).
bool Foam::primitiveMesh::checkFaceFlatness
(
+ const pointField& points,
+ const vectorField& faceCentres,
+ const vectorField& faceAreas,
const bool report,
const scalar warnFlatness,
labelHashSet* setPtr
@@ -992,52 +793,34 @@ bool Foam::primitiveMesh::checkFaceFlatness
<< exit(FatalError);
}
-
- const pointField& p = points();
const faceList& fcs = faces();
- const pointField& fctrs = faceCentres();
- // Areas are calculated as the sum of areas. (see
- // primitiveMeshFaceCentresAndAreas.C)
- scalarField magAreas(mag(faceAreas()));
+ tmp<scalarField> tfaceFlatness = primitiveMeshTools::faceFlatness
+ (
+ *this,
+ points,
+ faceCentres,
+ faceAreas
+ );
+ const scalarField& faceFlatness = tfaceFlatness();
- label nWarped = 0;
+ scalarField magAreas(mag(faceAreas));
scalar minFlatness = GREAT;
scalar sumFlatness = 0;
label nSummed = 0;
+ label nWarped = 0;
- forAll(fcs, faceI)
+ forAll(faceFlatness, faceI)
{
- const face& f = fcs[faceI];
-
- if (f.size() > 3 && magAreas[faceI] > VSMALL)
+ if (fcs[faceI].size() > 3 && magAreas[faceI] > VSMALL)
{
- 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);
- }
-
- scalar flatness = magAreas[faceI] / (sumA+VSMALL);
-
- sumFlatness += flatness;
+ sumFlatness += faceFlatness[faceI];
nSummed++;
- minFlatness = min(minFlatness, flatness);
+ minFlatness = min(minFlatness, faceFlatness[faceI]);
- if (flatness < warnFlatness)
+ if (faceFlatness[faceI] < warnFlatness)
{
nWarped++;
@@ -1092,131 +875,103 @@ bool Foam::primitiveMesh::checkFaceFlatness
}
-// Check 1D/2Dness of edges. Gets passed the non-empty directions and
-// checks all edges in the mesh whether they:
-// - have no component in a non-empty direction or
-// - are only in a singe non-empty direction.
-// Empty direction info is passed in as a vector of labels (synchronised)
-// which are 1 if the direction is non-empty, 0 if it is.
-bool Foam::primitiveMesh::checkEdgeAlignment
+bool Foam::primitiveMesh::checkConcaveCells
(
+ const vectorField& fAreas,
+ const pointField& fCentres,
const bool report,
- const Vector<label>& directions,
labelHashSet* setPtr
) const
{
if (debug)
{
- Info<< "bool primitiveMesh::checkEdgeAlignment("
- << "const bool, const Vector<label>&, labelHashSet*) const: "
- << "checking edge alignment" << endl;
+ Info<< "bool primitiveMesh::checkConcaveCells(const bool"
+ << ", labelHashSet*) const: "
+ << "checking for concave cells" << endl;
}
- label nDirs = 0;
- for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
- {
- if (directions[cmpt] == 1)
- {
- nDirs++;
- }
- else if (directions[cmpt] != 0)
- {
- FatalErrorIn
- (
- "primitiveMesh::checkEdgeAlignment"
- "(const bool, const Vector<label>&, labelHashSet*)"
- ) << "directions should contain 0 or 1 but is now " << directions
- << exit(FatalError);
- }
- }
+ const cellList& c = cells();
+ const labelList& fOwner = faceOwner();
- if (nDirs == vector::nComponents)
- {
- return false;
- }
+ label nConcaveCells = 0;
+ forAll(c, cellI)
+ {
+ const cell& cFaces = c[cellI];
+ bool concave = false;
- const pointField& p = points();
- const faceList& fcs = faces();
+ forAll(cFaces, i)
+ {
+ if (concave)
+ {
+ break;
+ }
- EdgeMap<label> edgesInError;
+ label fI = cFaces[i];
- forAll(fcs, faceI)
- {
- const face& f = fcs[faceI];
+ const point& fC = fCentres[fI];
- forAll(f, fp)
- {
- label p0 = f[fp];
- label p1 = f.nextLabel(fp);
- if (p0 < p1)
- {
- vector d(p[p1]-p[p0]);
- scalar magD = mag(d);
+ vector fN = fAreas[fI];
- if (magD > ROOTVSMALL)
- {
- d /= magD;
+ fN /= max(mag(fN), VSMALL);
- // Check how many empty directions are used by the edge.
- label nEmptyDirs = 0;
- label nNonEmptyDirs = 0;
- for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
- {
- if (mag(d[cmpt]) > 1e-6)
- {
- if (directions[cmpt] == 0)
- {
- nEmptyDirs++;
- }
- else
- {
- nNonEmptyDirs++;
- }
- }
- }
+ // Flip normal if required so that it is always pointing out of
+ // the cell
+ if (fOwner[fI] != cellI)
+ {
+ fN *= -1;
+ }
- if (nEmptyDirs == 0)
- {
- // Purely in ok directions.
- }
- else if (nEmptyDirs == 1)
+ // Is the centre of any other face of the cell on the
+ // wrong side of the plane of this face?
+
+ forAll(cFaces, j)
+ {
+ if (j != i)
+ {
+ label fJ = cFaces[j];
+
+ const point& pt = fCentres[fJ];
+
+ // If the cell is concave, the point will be on the
+ // positive normal side of the plane of f, defined by
+ // its centre and normal, and the angle between (pt -
+ // fC) and fN will be less than 90 degrees, so the dot
+ // product will be positive.
+
+ vector pC = (pt - fC);
+
+ pC /= max(mag(pC), VSMALL);
+
+ if ((pC & fN) > -planarCosAngle_)
{
- // Ok if purely in empty directions.
- if (nNonEmptyDirs > 0)
+ // Concave or planar face
+
+ concave = true;
+
+ if (setPtr)
{
- edgesInError.insert(edge(p0, p1), faceI);
+ setPtr->insert(cellI);
}
- }
- else if (nEmptyDirs > 1)
- {
- // Always an error
- edgesInError.insert(edge(p0, p1), faceI);
+
+ nConcaveCells++;
+
+ break;
}
}
}
}
}
- label nErrorEdges = returnReduce(edgesInError.size(), sumOp<label>());
+ reduce(nConcaveCells, sumOp<label>());
- if (nErrorEdges > 0)
+ if (nConcaveCells > 0)
{
if (debug || report)
{
- Info<< " ***Number of edges not aligned with or perpendicular to "
- << "non-empty directions: " << nErrorEdges << endl;
- }
-
- if (setPtr)
- {
- setPtr->resize(2*edgesInError.size());
- forAllConstIter(EdgeMap<label>, edgesInError, iter)
- {
- setPtr->insert(iter.key()[0]);
- setPtr->insert(iter.key()[1]);
- }
+ Info<< " ***Concave cells (using face planes) found,"
+ << " number of cells: " << nConcaveCells << endl;
}
return true;
@@ -1225,14 +980,18 @@ bool Foam::primitiveMesh::checkEdgeAlignment
{
if (debug || report)
{
- Info<< " All edges aligned with or perpendicular to "
- << "non-empty directions." << endl;
+ Info<< " Concave cell check OK." << endl;
}
+
return false;
}
+
+ return false;
}
+// Topological tests
+
bool Foam::primitiveMesh::checkUpperTriangular
(
const bool report,
@@ -1574,6 +1333,80 @@ bool Foam::primitiveMesh::checkFaceVertices
}
+bool Foam::primitiveMesh::checkPoints
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ if (debug)
+ {
+ Info<< "bool primitiveMesh::checkPoints"
+ << "(const bool, labelHashSet*) const: "
+ << "checking points" << endl;
+ }
+
+ label nFaceErrors = 0;
+ label nCellErrors = 0;
+
+ const labelListList& pf = pointFaces();
+
+ forAll(pf, pointI)
+ {
+ if (pf[pointI].empty())
+ {
+ if (setPtr)
+ {
+ setPtr->insert(pointI);
+ }
+
+ nFaceErrors++;
+ }
+ }
+
+
+ forAll(pf, pointI)
+ {
+ const labelList& pc = pointCells(pointI);
+
+ if (pc.empty())
+ {
+ if (setPtr)
+ {
+ setPtr->insert(pointI);
+ }
+
+ nCellErrors++;
+ }
+ }
+
+ reduce(nFaceErrors, sumOp<label>());
+ reduce(nCellErrors, sumOp<label>());
+
+ if (nFaceErrors > 0 || nCellErrors > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " ***Unused points found in the mesh, "
+ "number unused by faces: " << nFaceErrors
+ << " number unused by cells: " << nCellErrors
+ << endl;
+ }
+
+ return true;
+ }
+ else
+ {
+ if (debug || report)
+ {
+ Info<< " Point usage OK." << endl;
+ }
+
+ return false;
+ }
+}
+
+
// Check if all points on face are shared between faces.
bool Foam::primitiveMesh::checkDuplicateFaces
(
@@ -1891,168 +1724,153 @@ bool Foam::primitiveMesh::checkFaceFaces
}
-// Checks cells with 1 or less internal faces. Give numerical problems.
-bool Foam::primitiveMesh::checkCellDeterminant
-(
- const bool report, // report,
- labelHashSet* setPtr, // setPtr
- const Vector<label>& meshD
-) const
-{
- if (debug)
- {
- Info<< "bool primitiveMesh::checkCellDeterminant(const bool"
- << ", labelHashSet*) const: "
- << "checking for under-determined cells" << endl;
- }
-
- // 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;
- }
- }
-
-
- const cellList& c = cells();
-
- label nErrorCells = 0;
-
- scalar minDet = GREAT;
- scalar sumDet = 0;
- label nSummed = 0;
-
- if (nDims == 1)
- {
- minDet = 1;
- sumDet = c.size()*minDet;
- nSummed = c.size();
- }
- else
- {
- forAll (c, cellI)
- {
- const labelList& curFaces = c[cellI];
-
- // Calculate local normalization factor
- scalar avgArea = 0;
-
- label nInternalFaces = 0;
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
- forAll(curFaces, i)
- {
- if (isInternalFace(curFaces[i]))
- {
- avgArea += mag(faceAreas()[curFaces[i]]);
+bool Foam::primitiveMesh::checkClosedBoundary(const bool report) const
+{
+ return checkClosedBoundary(faceAreas(), report, PackedBoolList(0));
+}
- nInternalFaces++;
- }
- }
- if (nInternalFaces == 0)
- {
- if (setPtr)
- {
- setPtr->insert(cellI);
- }
+bool Foam::primitiveMesh::checkClosedCells
+(
+ const bool report,
+ labelHashSet* setPtr,
+ labelHashSet* aspectSetPtr,
+ const Vector<label>& solutionD
+) const
+{
+ return checkClosedCells
+ (
+ faceAreas(),
+ cellVolumes(),
+ report,
+ setPtr,
+ aspectSetPtr,
+ solutionD
+ );
+}
- nErrorCells++;
- }
- else
- {
- avgArea /= nInternalFaces;
- symmTensor areaTensor(symmTensor::zero);
+bool Foam::primitiveMesh::checkFaceAreas
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ return checkFaceAreas
+ (
+ faceAreas(),
+ report,
+ false, // detailedReport,
+ setPtr
+ );
+}
- forAll(curFaces, i)
- {
- if (isInternalFace(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;
- }
- }
+bool Foam::primitiveMesh::checkCellVolumes
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ return checkCellVolumes
+ (
+ cellVolumes(),
+ report,
+ false, // detailedReport,
+ setPtr
+ );
+}
- scalar determinant = mag(det(areaTensor));
- minDet = min(determinant, minDet);
- sumDet += determinant;
- nSummed++;
+bool Foam::primitiveMesh::checkFaceOrthogonality
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ return checkFaceOrthogonality
+ (
+ faceAreas(),
+ cellCentres(),
+ report,
+ setPtr
+ );
+}
- if (determinant < 1e-3)
- {
- if (setPtr)
- {
- setPtr->insert(cellI);
- }
- nErrorCells++;
- }
- }
- }
- }
+bool Foam::primitiveMesh::checkFacePyramids
+(
+ const bool report,
+ const scalar minPyrVol,
+ labelHashSet* setPtr
+) const
+{
+ return checkFacePyramids
+ (
+ points(),
+ cellCentres(),
+ report,
+ false, // detailedReport,
+ minPyrVol,
+ setPtr
+ );
+}
- reduce(nErrorCells, sumOp<label>());
- reduce(minDet, minOp<scalar>());
- reduce(sumDet, sumOp<scalar>());
- reduce(nSummed, sumOp<label>());
- if (debug || report)
- {
- if (nSummed > 0)
- {
- Info<< " Cell determinant (wellposedness) : minimum: " << minDet
- << " average: " << sumDet/nSummed
- << endl;
- }
- }
+bool Foam::primitiveMesh::checkFaceSkewness
+(
+ const bool report,
+ labelHashSet* setPtr
+) const
+{
+ return checkFaceSkewness
+ (
+ points(),
+ faceCentres(),
+ faceAreas(),
+ cellCentres(),
+ report,
+ setPtr
+ );
+}
- if (nErrorCells > 0)
- {
- if (debug || report)
- {
- Info<< " ***Cells with small determinant found, number of cells: "
- << nErrorCells << endl;
- }
- return true;
- }
- else
- {
- if (debug || report)
- {
- Info<< " Cell determinant check OK." << endl;
- }
+bool Foam::primitiveMesh::checkFaceAngles
+(
+ const bool report,
+ const scalar maxDeg,
+ labelHashSet* setPtr
+) const
+{
+ return checkFaceAngles
+ (
+ points(),
+ faceAreas(),
+ report,
+ maxDeg,
+ setPtr
+ );
+}
- return false;
- }
- return false;
+bool Foam::primitiveMesh::checkFaceFlatness
+(
+ const bool report,
+ const scalar warnFlatness,
+ labelHashSet* setPtr
+) const
+{
+ return checkFaceFlatness
+ (
+ points(),
+ faceCentres(),
+ faceAreas(),
+ report,
+ warnFlatness,
+ setPtr
+ );
}
@@ -2062,117 +1880,16 @@ bool Foam::primitiveMesh::checkConcaveCells
labelHashSet* setPtr
) const
{
- if (debug)
- {
- Info<< "bool primitiveMesh::checkConcaveCells(const bool"
- << ", labelHashSet*) const: "
- << "checking for concave cells" << endl;
- }
-
- const cellList& c = cells();
- const labelList& fOwner = faceOwner();
- const vectorField& fAreas = faceAreas();
- const pointField& fCentres = faceCentres();
-
- label nConcaveCells = 0;
-
- forAll(c, cellI)
- {
- const cell& cFaces = c[cellI];
-
- bool concave = false;
-
- forAll(cFaces, i)
- {
- if (concave)
- {
- break;
- }
-
- label fI = cFaces[i];
-
- const point& fC = fCentres[fI];
-
- vector fN = fAreas[fI];
-
- fN /= max(mag(fN), VSMALL);
-
- // Flip normal if required so that it is always pointing out of
- // the cell
- if (fOwner[fI] != cellI)
- {
- fN *= -1;
- }
-
- // Is the centre of any other face of the cell on the
- // wrong side of the plane of this face?
-
- forAll(cFaces, j)
- {
- if (j != i)
- {
- label fJ = cFaces[j];
-
- const point& pt = fCentres[fJ];
-
- // If the cell is concave, the point will be on the
- // positive normal side of the plane of f, defined by
- // its centre and normal, and the angle between (pt -
- // fC) and fN will be less than 90 degrees, so the dot
- // product will be positive.
-
- vector pC = (pt - fC);
-
- pC /= max(mag(pC), VSMALL);
-
- if ((pC & fN) > -planarCosAngle_)
- {
- // Concave or planar face
-
- concave = true;
-
- if (setPtr)
- {
- setPtr->insert(cellI);
- }
-
- nConcaveCells++;
-
- break;
- }
- }
- }
- }
- }
-
- reduce(nConcaveCells, sumOp<label>());
-
- if (nConcaveCells > 0)
- {
- if (debug || report)
- {
- Info<< " ***Concave cells (using face planes) found,"
- << " number of cells: " << nConcaveCells << endl;
- }
-
- return true;
- }
- else
- {
- if (debug || report)
- {
- Info<< " Concave cell check OK." << endl;
- }
-
- return false;
- }
-
- return false;
+ return checkConcaveCells
+ (
+ faceAreas(),
+ faceCentres(),
+ report,
+ setPtr
+ );
}
-// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-
bool Foam::primitiveMesh::checkTopology(const bool report) const
{
label noFailedChecks = 0;
@@ -2182,7 +1899,6 @@ bool Foam::primitiveMesh::checkTopology(const bool report) const
if (checkCellsZipUp(report)) noFailedChecks++;
if (checkFaceVertices(report)) noFailedChecks++;
if (checkFaceFaces(report)) noFailedChecks++;
- //if (checkCellDeterminant(report)) noFailedChecks++;
if (noFailedChecks == 0)
{
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C
deleted file mode 100644
index 4b0d2da63aa25fb48475dca70dca108eea3d0354..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C
+++ /dev/null
@@ -1,266 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration |
- \\ / A nd | Copyright (C) 2011 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
- Given a set of points, find out if the mesh resulting from point motion will
- be valid without actually changing the mesh.
-
-\*---------------------------------------------------------------------------*/
-
-#include "primitiveMesh.H"
-#include "pyramidPointFaceRef.H"
-#include "cell.H"
-#include "unitConversion.H"
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-bool Foam::primitiveMesh::checkMeshMotion
-(
- const pointField& newPoints,
- const bool report
-) const
-{
- if (debug || report)
- {
- Pout<< "bool primitiveMesh::checkMeshMotion("
- << "const pointField& newPoints, const bool report) const: "
- << "checking mesh motion" << endl;
- }
-
- bool error = false;
-
- const faceList& f = faces();
-
- const labelList& own = faceOwner();
- const labelList& nei = faceNeighbour();
-
- vectorField fCtrs(nFaces());
- vectorField fAreas(nFaces());
-
- makeFaceCentresAndAreas(newPoints, fCtrs, fAreas);
-
- // Check cell volumes and calculate new cell centres
- vectorField cellCtrs(nCells());
- scalarField cellVols(nCells());
-
- makeCellCentresAndVols(fCtrs, fAreas, cellCtrs, cellVols);
-
- scalar minVolume = GREAT;
- label nNegVols = 0;
-
- forAll(cellVols, cellI)
- {
- if (cellVols[cellI] < VSMALL)
- {
- if (debug || report)
- {
- Pout<< "Zero or negative cell volume detected for cell "
- << cellI << ". Volume = " << cellVols[cellI] << endl;
- }
-
- nNegVols++;
- }
-
- minVolume = min(minVolume, cellVols[cellI]);
- }
-
- if (nNegVols > 0)
- {
- error = true;
-
- Pout<< "Zero or negative cell volume in mesh motion in " << nNegVols
- << " cells. Min volume: " << minVolume << endl;
- }
- else
- {
- if (debug || report)
- {
- Pout<< "Min volume = " << minVolume
- << ". Total volume = " << sum(cellVols)
- << ". Cell volumes OK." << endl;
- }
- }
-
- // Check face areas
-
- scalar minArea = GREAT;
- label nNegAreas = 0;
- label nPyrErrors = 0;
- label nDotProductErrors = 0;
-
- forAll(f, faceI)
- {
- const scalar a = Foam::mag(fAreas[faceI]);
-
- if (a < VSMALL)
- {
- if (debug || report)
- {
- if (isInternalFace(faceI))
- {
- Pout<< "Zero or negative face area detected for "
- << "internal face "<< faceI << " between cells "
- << own[faceI] << " and " << nei[faceI]
- << ". Face area magnitude = " << a << endl;
- }
- else
- {
- Pout<< "Zero or negative face area detected for "
- << "boundary face " << faceI << " next to cell "
- << own[faceI] << ". Face area magnitude = "
- << a << endl;
- }
- }
-
- nNegAreas++;
- }
-
- minArea = min(minArea, a);
-
- // Create the owner pyramid - it will have negative volume
- scalar pyrVol =
- pyramidPointFaceRef(f[faceI], cellCtrs[own[faceI]]).mag(newPoints);
-
- if (pyrVol > SMALL)
- {
- if (debug || report)
- {
- Pout<< "Negative pyramid volume: " << -pyrVol
- << " for face " << faceI << " " << f[faceI]
- << " and owner cell: " << own[faceI] << endl
- << "Owner cell vertex labels: "
- << cells()[own[faceI]].labels(f)
- << endl;
- }
-
- nPyrErrors++;
- }
-
- if (isInternalFace(faceI))
- {
- // Create the neighbour pyramid - it will have positive volume
- scalar pyrVol =
- pyramidPointFaceRef
- (
- f[faceI],
- cellCtrs[nei[faceI]]
- ).mag(newPoints);
-
- if (pyrVol < -SMALL)
- {
- if (debug || report)
- {
- Pout<< "Negative pyramid volume: " << pyrVol
- << " for face " << faceI << " " << f[faceI]
- << " and neighbour cell: " << nei[faceI] << nl
- << "Neighbour cell vertex labels: "
- << cells()[nei[faceI]].labels(f)
- << endl;
- }
-
- nPyrErrors++;
- }
-
- const vector d = cellCtrs[nei[faceI]] - cellCtrs[own[faceI]];
- const vector& s = fAreas[faceI];
- scalar dDotS = (d & s)/(mag(d)*mag(s) + VSMALL);
-
- // Only write full message the first time
- if (dDotS < SMALL && nDotProductErrors == 0)
- {
- // Non-orthogonality greater than 90 deg
- WarningIn
- (
- "primitiveMesh::checkMeshMotion"
- "(const pointField& newPoints, const bool report) const"
- ) << "Severe non-orthogonality in mesh motion for face "
- << faceI
- << " between cells " << own[faceI] << " and " << nei[faceI]
- << ": Angle = " << radToDeg(::acos(dDotS))
- << " deg." << endl;
-
- nDotProductErrors++;
- }
- }
- }
-
- if (nNegAreas > 0)
- {
- error = true;
-
- WarningIn
- (
- "primitiveMesh::checkMeshMotion"
- "(const pointField& newPoints, const bool report) const"
- ) << "Zero or negative face area in mesh motion in " << nNegAreas
- << " faces. Min area: " << minArea << endl;
- }
- else
- {
- if (debug || report)
- {
- Pout<< "Min area = " << minArea
- << ". Face areas OK." << endl;
- }
- }
-
- if (nPyrErrors > 0)
- {
- Pout<< "Detected " << nPyrErrors
- << " negative pyramid volume in mesh motion" << endl;
-
- error = true;
- }
- else
- {
- if (debug || report)
- {
- Pout<< "Pyramid volumes OK." << endl;
- }
- }
-
- if (nDotProductErrors > 0)
- {
- Pout<< "Detected " << nDotProductErrors
- << " in non-orthogonality in mesh motion." << endl;
-
- error = true;
- }
- else
- {
- if (debug || report)
- {
- Pout<< "Non-orthogonality check OK." << endl;
- }
- }
-
- if (!error && (debug || report))
- {
- Pout<< "Mesh motion check OK." << endl;
- }
-
- return error;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
new file mode 100644
index 0000000000000000000000000000000000000000..e94c8edbf77a1709d63377a11c6fb98e5765e067
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
@@ -0,0 +1,638 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 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::primitiveMeshTools::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) + SMALL))*d;
+ vector svHat = sv/(mag(sv) + VSMALL);
+
+ // 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) + VSMALL;
+ 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::primitiveMeshTools::faceOrthogonality
+(
+ const point& ownCc,
+ const point& neiCc,
+ const vector& s
+)
+{
+ vector d = neiCc - ownCc;
+
+ return (d & s)/(mag(d)*mag(s) + VSMALL);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::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();
+
+ // Internal faces
+ forAll(nei, faceI)
+ {
+ ortho[faceI] = faceOrthogonality
+ (
+ cc[own[faceI]],
+ cc[nei[faceI]],
+ areas[faceI]
+ );
+ }
+
+ return tortho;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::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();
+ const faceList& fcs = mesh.faces();
+
+ tmp<scalarField> tskew(new scalarField(mesh.nFaces()));
+ scalarField& skew = tskew();
+
+ 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++)
+ {
+ vector Cpf = fCtrs[faceI] - cellCtrs[own[faceI]];
+
+ vector normal = fAreas[faceI];
+ normal /= mag(normal) + VSMALL;
+ vector d = normal*(normal & Cpf);
+
+
+ // Skewness vector
+ vector sv =
+ Cpf
+ - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
+ vector svHat = sv/(mag(sv) + VSMALL);
+
+ // 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) + VSMALL;
+ const face& f = fcs[faceI];
+ forAll(f, pi)
+ {
+ fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
+ }
+
+ // Normalised skewness
+ skew[faceI] = mag(sv)/fd;
+ }
+
+ return tskew;
+}
+
+void Foam::primitiveMeshTools::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::primitiveMeshTools::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(), vector::zero);
+ vectorField sumMagClosed(mesh.nCells(), vector::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] + VSMALL)
+ );
+ }
+ 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 + VSMALL);
+ if (nDims == 3)
+ {
+ aspectRatio = max
+ (
+ aspectRatio,
+ 1.0/6.0*cmptSum(sumMagClosed[cellI])/pow(vols[cellI], 2.0/3.0)
+ );
+ }
+
+ aratio[cellI] = aspectRatio;
+ }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceConcavity
+(
+ const scalar maxSin,
+ const primitiveMesh& mesh,
+ const pointField& p,
+ const vectorField& faceAreas
+)
+{
+ const faceList& fcs = mesh.faces();
+
+ vectorField faceNormals(faceAreas);
+ faceNormals /= mag(faceNormals) + VSMALL;
+
+ tmp<scalarField> tfaceAngles(new scalarField(mesh.nFaces()));
+ scalarField& faceAngles = tfaceAngles();
+
+
+ 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 + VSMALL;
+
+ scalar maxEdgeSin = 0.0;
+
+ forAll(f, fp0)
+ {
+ // Get vertex after fp
+ label fp1 = f.fcIndex(fp0);
+
+ // Normalized vector between two consecutive points
+ vector e10(p[f[fp1]] - p[f[fp0]]);
+ scalar magE10 = mag(e10);
+ e10 /= magE10 + VSMALL;
+
+ 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::primitiveMeshTools::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();
+
+
+ forAll(fcs, faceI)
+ {
+ const face& f = fcs[faceI];
+
+ if (f.size() > 3 && magAreas[faceI] > VSMALL)
+ {
+ 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+VSMALL);
+ }
+ }
+
+ return tfaceFlatness;
+}
+
+
+//Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::edgeAlignment
+//(
+// const primitiveMesh& mesh,
+// const Vector<label>& meshD,
+// const pointField& p
+//)
+//{
+// label nDirs = 0;
+// for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
+// {
+// if (meshD[cmpt] == 1)
+// {
+// nDirs++;
+// }
+// else if (meshD[cmpt] != 0)
+// {
+// FatalErrorIn
+// (
+// "primitiveMeshTools::edgeAlignment"
+// "(const primitiveMesh&, const Vector<label>&"
+// ", const pointField&)"
+// ) << "directions should contain 0 or 1 but is now " << meshD
+// << exit(FatalError);
+// }
+// }
+//
+// tmp<scalarField> tedgeAlignment(new scalarField(mesh.nFaces(), 0.0));
+// scalarField& edgeAlignment = tedgeAlignment();
+//
+// if (nDirs != vector::nComponents)
+// {
+// const faceList& fcs = mesh.faces();
+//
+// forAll(fcs, faceI)
+// {
+// const face& f = fcs[faceI];
+//
+// forAll(f, fp)
+// {
+// label p0 = f[fp];
+// label p1 = f.nextLabel(fp);
+// if (p0 < p1)
+// {
+// vector d(p[p1]-p[p0]);
+// scalar magD = mag(d);
+//
+// if (magD > ROOTVSMALL)
+// {
+// d /= magD;
+//
+// // Check how many empty directions are used by the
+// // edge.
+// label nEmptyDirs = 0;
+// label nNonEmptyDirs = 0;
+// for
+// (
+// direction cmpt=0;
+// cmpt<vector::nComponents;
+// cmpt++
+// )
+// {
+// if (mag(d[cmpt]) > 1e-6)
+// {
+// if (meshD[cmpt] == 0)
+// {
+// nEmptyDirs++;
+// }
+// else
+// {
+// nNonEmptyDirs++;
+// }
+// }
+// }
+//
+// if (nEmptyDirs == 0)
+// {
+// // Purely in ok directions.
+// }
+// else if (nEmptyDirs == 1)
+// {
+// // Ok if purely in empty directions.
+// if (nNonEmptyDirs > 0)
+// {
+// edgeAlignment[faceI] = max
+// (
+// edgeAlignment[faceI],
+// 1
+// );
+// }
+// }
+// else if (nEmptyDirs > 1)
+// {
+// // Always an error
+// edgeAlignment[faceI] = max
+// (
+// edgeAlignment[faceI],
+// 2
+// );
+// }
+// }
+// }
+// }
+// }
+// }
+//
+// return tedgeAlignment;
+//}
+
+
+// Checks cells with 1 or less internal faces. Give numerical problems.
+Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::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();
+
+ const cellList& c = mesh.cells();
+
+ if (nDims == 1)
+ {
+ cellDeterminant = 1.0;
+ }
+ else
+ {
+ forAll (c, cellI)
+ {
+ const labelList& curFaces = c[cellI];
+
+ // Calculate local normalization 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(symmTensor::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;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H
new file mode 100644
index 0000000000000000000000000000000000000000..5e3b2c07fd4b03aaacab2f82fca09c83187e74c9
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H
@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2012 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::primitiveMeshTools
+
+Description
+ Collection of static functions operating on primitiveMesh (mainly checks).
+
+SourceFiles
+ primitiveMeshTools.C
+
+\*---------------------------------------------------------------------------*/
+#ifndef primitiveMeshTools_H
+#define primitiveMeshTools_H
+
+#include "primitiveMesh.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Namespace primitiveMeshTools Declaration
+\*---------------------------------------------------------------------------*/
+
+class primitiveMeshTools
+{
+
+protected:
+
+ // Protected Member Functions
+
+ //- Skewness of single face
+ static scalar faceSkewness
+ (
+ const primitiveMesh& mesh,
+ const pointField& p,
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+
+ const label faceI,
+ const point& ownCc,
+ const point& neiCc
+ );
+
+ //- Orthogonality of single face
+ static scalar faceOrthogonality
+ (
+ const point& ownCc,
+ const point& neiCc,
+ const vector& s
+ );
+
+public:
+
+ //- Generate non-orthogonality field (internal faces only)
+ static tmp<scalarField> faceOrthogonality
+ (
+ const primitiveMesh& mesh,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs
+ );
+
+ //- Generate face pyramid volume fields
+ static void facePyramidVolume
+ (
+ const primitiveMesh& mesh,
+ const pointField& points,
+ const vectorField& cellCtrs,
+ scalarField& ownPyrVol,
+ scalarField& neiPyrVol
+ );
+
+ //- Generate skewness field
+ static tmp<scalarField> faceSkewness
+ (
+ const primitiveMesh& mesh,
+ const pointField& points,
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs
+ );
+
+ //- Generate cell openness and cell ascpect ratio field
+ static 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
+ static 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)
+ static 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)
+ static tmp<scalarField> edgeAlignment
+ (
+ const primitiveMesh& mesh,
+ const Vector<label>& directions,
+ const pointField& p
+ );
+
+ //- Generate cell determinant field
+ static tmp<scalarField> cellDeterminant
+ (
+ const primitiveMesh& mesh,
+ const Vector<label>& directions,
+ const vectorField& faceAreas,
+ const PackedBoolList& internalOrCoupledFace
+ );
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //