From e58ff9d366c4f70aa044a0a199ef71a0ceee5400 Mon Sep 17 00:00:00 2001
From: Henry <Henry>
Date: Thu, 5 Feb 2015 16:31:29 +0000
Subject: [PATCH] primitiveMeshTools: stabilize with VSMALL rather than SMALL
 to avoid problems with very small meshes Resolves bug-report
 http://www.openfoam.org/mantisbt/view.php?id=1509

---
 .../primitiveMeshCheck/primitiveMeshTools.C   | 120 +-----------------
 1 file changed, 5 insertions(+), 115 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
index fa83d3ae9d..868c53954c 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,7 +47,7 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness
     // Skewness vector
     vector sv =
         Cpf
-      - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + SMALL))*d;
+      - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + VSMALL))*d;
     vector svHat = sv/(mag(sv) + VSMALL);
 
     // Normalisation distance calculated as the approximate distance
@@ -64,6 +64,7 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness
     return mag(sv)/fd;
 }
 
+
 Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness
 (
     const primitiveMesh& mesh,
@@ -196,6 +197,7 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceSkewness
     return tskew;
 }
 
+
 void Foam::primitiveMeshTools::facePyramidVolume
 (
     const primitiveMesh& mesh,
@@ -440,7 +442,7 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceFlatness
                 sumA += mag(n);
             }
 
-            faceFlatness[faceI] = magAreas[faceI] / (sumA+VSMALL);
+            faceFlatness[faceI] = magAreas[faceI]/(sumA + VSMALL);
         }
     }
 
@@ -448,116 +450,6 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::faceFlatness
 }
 
 
-//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,
@@ -581,8 +473,6 @@ Foam::tmp<Foam::scalarField> Foam::primitiveMeshTools::cellDeterminant
         }
     }
 
-
-
     tmp<scalarField> tcellDeterminant(new scalarField(mesh.nCells()));
     scalarField& cellDeterminant = tcellDeterminant();