diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 76c5f1c8a3777cc65b520d2a720fc1fc66ddd895..a416fafc301d52778034467bdf9f83e649956d08 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
@@ -265,6 +265,24 @@ private:
const Vector<label>& meshD
) const;
+ bool checkFaceWeight
+ (
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
+ const bool report,
+ const scalar minWeight,
+ labelHashSet* setPtr
+ ) const;
+
+ bool checkVolRatio
+ (
+ const scalarField& cellVols,
+ const bool report,
+ const scalar minRatio,
+ labelHashSet* setPtr
+ ) const;
+
public:
// Public typedefs
@@ -612,6 +630,22 @@ public:
const bool detailedReport = false
) const;
+ //- Check for face weights
+ virtual bool checkFaceWeight
+ (
+ const bool report,
+ const scalar minWeight = 0.05,
+ labelHashSet* setPtr = NULL
+ ) const;
+
+ //- Check for neighbouring cell volumes
+ virtual bool checkVolRatio
+ (
+ const bool report,
+ const scalar minRatio = 0.01,
+ labelHashSet* setPtr = NULL
+ ) const;
+
// Helper functions
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
index f7daba04955b96710f570267d344519909a3d9f9..45bb3d67e475a3efe3f6ca6d8404c4cea7cf357a 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshCheck.C
@@ -500,6 +500,190 @@ bool Foam::polyMesh::checkCellDeterminant
}
+bool Foam::polyMesh::checkFaceWeight
+(
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs,
+ const bool report,
+ const scalar minWeight,
+ labelHashSet* setPtr
+) const
+{
+ if (debug)
+ {
+ Info<< "bool polyMesh::checkFaceWeight(const bool"
+ << ", labelHashSet*) const: "
+ << "checking for low face interpolation weights" << endl;
+ }
+
+ tmp<scalarField> tfaceWght = polyMeshTools::faceWeights
+ (
+ *this,
+ fCtrs,
+ fAreas,
+ cellCtrs
+ );
+ scalarField& faceWght = tfaceWght();
+
+
+ label nErrorFaces = 0;
+ scalar minDet = GREAT;
+ scalar sumDet = 0.0;
+ label nSummed = 0;
+
+ // Statistics only for internal and masters of coupled faces
+ PackedBoolList isMasterFace(syncTools::getInternalOrMasterFaces(*this));
+
+ forAll (faceWght, faceI)
+ {
+ if (faceWght[faceI] < minWeight)
+ {
+ // Note: insert both sides of coupled faces
+ if (setPtr)
+ {
+ setPtr->insert(faceI);
+ }
+
+ nErrorFaces++;
+ }
+
+ // Note: statistics only on master of coupled faces
+ if (isMasterFace[faceI])
+ {
+ minDet = min(minDet, faceWght[faceI]);
+ sumDet += faceWght[faceI];
+ nSummed++;
+ }
+ }
+
+ reduce(nErrorFaces, sumOp<label>());
+ reduce(minDet, minOp<scalar>());
+ reduce(sumDet, sumOp<scalar>());
+ reduce(nSummed, sumOp<label>());
+
+ if (debug || report)
+ {
+ if (nSummed > 0)
+ {
+ Info<< " Face interpolation weight : minimum: " << minDet
+ << " average: " << sumDet/nSummed
+ << endl;
+ }
+ }
+
+ if (nErrorFaces > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " ***Faces with small interpolation found, number of faces: "
+ << nErrorFaces << endl;
+ }
+
+ return true;
+ }
+ else
+ {
+ if (debug || report)
+ {
+ Info<< " Face interpolation weight check OK." << endl;
+ }
+
+ return false;
+ }
+
+ return false;
+}
+
+
+bool Foam::polyMesh::checkVolRatio
+(
+ const scalarField& cellVols,
+ const bool report,
+ const scalar minRatio,
+ labelHashSet* setPtr
+) const
+{
+ if (debug)
+ {
+ Info<< "bool polyMesh::checkVolRatio(const bool"
+ << ", labelHashSet*) const: "
+ << "checking for low volume ratio" << endl;
+ }
+
+ tmp<scalarField> tvolRatio = polyMeshTools::volRatio(*this, cellVols);
+ scalarField& volRatio = tvolRatio();
+
+
+ label nErrorFaces = 0;
+ scalar minDet = GREAT;
+ scalar sumDet = 0.0;
+ label nSummed = 0;
+
+ // Statistics only for internal and masters of coupled faces
+ PackedBoolList isMasterFace(syncTools::getInternalOrMasterFaces(*this));
+
+ forAll (volRatio, faceI)
+ {
+ if (volRatio[faceI] < minRatio)
+ {
+ // Note: insert both sides of coupled faces
+ if (setPtr)
+ {
+ setPtr->insert(faceI);
+ }
+
+ nErrorFaces++;
+ }
+
+ // Note: statistics only on master of coupled faces
+ if (isMasterFace[faceI])
+ {
+ minDet = min(minDet, volRatio[faceI]);
+ sumDet += volRatio[faceI];
+ nSummed++;
+ }
+ }
+
+ reduce(nErrorFaces, sumOp<label>());
+ reduce(minDet, minOp<scalar>());
+ reduce(sumDet, sumOp<scalar>());
+ reduce(nSummed, sumOp<label>());
+
+ if (debug || report)
+ {
+ if (nSummed > 0)
+ {
+ Info<< " Face volume ratio : minimum: " << minDet
+ << " average: " << sumDet/nSummed
+ << endl;
+ }
+ }
+
+ if (nErrorFaces > 0)
+ {
+ if (debug || report)
+ {
+ Info<< " ***Faces with small volume ratio found, number of faces: "
+ << nErrorFaces << endl;
+ }
+
+ return true;
+ }
+ else
+ {
+ if (debug || report)
+ {
+ Info<< " Face volume ratio check OK." << endl;
+ }
+
+ return false;
+ }
+
+ return false;
+}
+
+
//- Could override checkClosedBoundary to not look at (collocated!) coupled
// faces
//bool Foam::polyMesh::checkClosedBoundary(const bool report) const
@@ -582,6 +766,36 @@ bool Foam::polyMesh::checkCellDeterminant
}
+bool Foam::polyMesh::checkFaceWeight
+(
+ const bool report,
+ const scalar minWeight,
+ labelHashSet* setPtr
+) const
+{
+ return checkFaceWeight
+ (
+ faceCentres(),
+ faceAreas(),
+ cellCentres(),
+ report,
+ minWeight,
+ setPtr
+ );
+}
+
+
+bool Foam::polyMesh::checkVolRatio
+(
+ const bool report,
+ const scalar minRatio,
+ labelHashSet* setPtr
+) const
+{
+ return checkVolRatio(cellVolumes(), report, minRatio, setPtr);
+}
+
+
bool Foam::polyMesh::checkMeshMotion
(
const pointField& newPoints,
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
index 55f6a907676e7f241fd34c11befaa23a5de4332f..1884485a69774f1916b80ee255f84dd440666999 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C
@@ -187,4 +187,112 @@ Foam::tmp<Foam::scalarField> Foam::polyMeshTools::faceSkewness
}
+Foam::tmp<Foam::scalarField> Foam::polyMeshTools::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();
+
+ // 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::polyMeshTools::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();
+
+ // Internal faces
+ forAll(nei, faceI)
+ {
+ scalar volOwn = vol[own[faceI]];
+ scalar volNei = vol[nei[faceI]];
+
+ ratio[faceI] = min(volOwn,volNei)/(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)/(volOwn+volNei+VSMALL);
+ }
+ }
+ }
+
+ return tratio;
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
index 53c758c3be7967f5a914b39e7e0c12a50bbaf86c..7f23facd04ec02eb0b43abb1b64ee2702f5be88e 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2013 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -95,6 +95,22 @@ public:
// );
// }
+ //- Generate interpolation factors field
+ static tmp<scalarField> faceWeights
+ (
+ const polyMesh& mesh,
+ const vectorField& fCtrs,
+ const vectorField& fAreas,
+ const vectorField& cellCtrs
+ );
+
+ //- Generate volume ratio field
+ static tmp<scalarField> volRatio
+ (
+ const polyMesh& mesh,
+ const scalarField& vol
+ );
+
};