From 8460d967cc64910229a98cd8a723f6827979bffb Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 15 Apr 2013 15:24:23 +0100
Subject: [PATCH 1/2] ENH: polyMesh: additional mesh checks about interpolation
 weights and volume ratio

---
 src/OpenFOAM/meshes/polyMesh/polyMesh.H       |  34 +++
 .../polyMesh/polyMeshCheck/polyMeshCheck.C    | 214 ++++++++++++++++++
 .../polyMesh/polyMeshCheck/polyMeshTools.C    | 108 +++++++++
 .../polyMesh/polyMeshCheck/polyMeshTools.H    |  18 +-
 4 files changed, 373 insertions(+), 1 deletion(-)

diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.H b/src/OpenFOAM/meshes/polyMesh/polyMesh.H
index 76c5f1c8a3..a416fafc30 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 f7daba0495..45bb3d67e4 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 55f6a90767..1884485a69 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 53c758c3be..7f23facd04 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
+    );
+
 };
 
 

From 415cdb6a6300623c7b847e4d9eda862125123406 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Mon, 15 Apr 2013 16:55:59 +0100
Subject: [PATCH 2/2] ENH: checkMesh: check volume ratio and face interpolation
 weights

---
 .../manipulation/checkMesh/checkGeometry.C    | 34 +++++++++++++++++++
 1 file changed, 34 insertions(+)

diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index ef521a257d..cf503b0574 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -832,5 +832,39 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
+    if (allGeometry)
+    {
+        faceSet faces(mesh, "lowWeightFaces", mesh.nFaces()/100);
+        if (mesh.checkFaceWeight(true, 0.05, &faces))
+        {
+            noFailedChecks++;
+
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nFaces
+                << " faces with low interpolation weights to set "
+                << faces.name() << endl;
+            faces.instance() = mesh.pointsInstance();
+            faces.write();
+        }
+    }
+
+    if (allGeometry)
+    {
+        faceSet faces(mesh, "lowVolRatioFaces", mesh.nFaces()/100);
+        if (mesh.checkVolRatio(true, 0.05, &faces))
+        {
+            noFailedChecks++;
+
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nFaces
+                << " faces with low volume ratio cells to set "
+                << faces.name() << endl;
+            faces.instance() = mesh.pointsInstance();
+            faces.write();
+        }
+    }
+
     return noFailedChecks;
 }