ENH: polyMesh: additional mesh checks about interpolation weights and volume ratio

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
 

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,

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;
+}
+
+
 // ************************************************************************* //

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
+    );
+
 };