From e1938c7a0cae8e2f0f4df5881f0a57b34df023b2 Mon Sep 17 00:00:00 2001
From: mattijs <m.janssens@opencfd.co.uk>
Date: Thu, 8 May 2008 23:42:20 +0100
Subject: [PATCH] Edge alignment check for 2D cases

---
 .../manipulation/checkMesh/checkGeometry.C    |  34 ++++-
 .../manipulation/checkMesh/checkTopology.C    |   4 +-
 .../meshes/primitiveMesh/primitiveMesh.H      |   9 ++
 .../primitiveMeshCheck/primitiveMeshCheck.C   | 141 ++++++++++++++++++
 4 files changed, 181 insertions(+), 7 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index c7b29edda40..6b1a10d90e6 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -4,6 +4,7 @@
 #include "cellSet.H"
 #include "faceSet.H"
 #include "pointSet.H"
+#include "EdgeMap.H"
 
 Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
 {
@@ -21,12 +22,33 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
     // Min length
     scalar minDistSqr = magSqr(1e-6*(globalBb.max() - globalBb.min()));
 
-    Info<< "    Mesh (non-empty) directions "
-        << (mesh.directions() + Vector<label>::one)/2 << endl;
-    
+    // Non-empty directions
+    const Vector<label> validDirs = (mesh.directions() + Vector<label>::one)/2;
+
+    Info<< "    Mesh (non-empty) directions " << validDirs << endl;
+
+    scalar nGeomDims = mesh.nGeometricD();
+
     Info<< "    Mesh (non-empty, non-wedge) dimensions "
-        << mesh.nGeometricD() << endl;
-    
+        << nGeomDims << endl;
+
+    if (nGeomDims < 3)
+    {
+        pointSet nonAlignedPoints(mesh, "nonAlignedEdges", mesh.nPoints()/100);
+
+        if (mesh.checkEdgeAlignment(true, validDirs, &nonAlignedPoints))
+        {
+            noFailedChecks++;
+
+            if (nonAlignedPoints.size() > 0)
+            {
+                Pout<< "  <<Writing " << nonAlignedPoints.size()
+                    << " points on non-aligned edges to set "
+                    << nonAlignedPoints.name() << endl;
+                nonAlignedPoints.write();
+            }
+        }
+    }
 
     if (mesh.checkClosedBoundary(true)) noFailedChecks++;
 
@@ -209,7 +231,7 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
             cells.write();
         }
     }
-    
+
 
     return noFailedChecks;
 }
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
index 2cb38266278..6fbb2c0bf9e 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkTopology.C
@@ -36,7 +36,9 @@ Foam::label Foam::checkTopology(const polyMesh& mesh, const bool allTopology)
         if (mesh.checkUpperTriangular(true, &faces))
         {
             noFailedChecks++;
-
+        }
+        if (faces.size() > 0)
+        {
             Pout<< "  <<Writing " << faces.size()
                 << " unordered faces to set " << faces.name() << endl;
             faces.write();
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
index 1dfc4c2e15c..18f0ddef9bd 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H
@@ -65,6 +65,7 @@ SourceFiles
 #include "boolList.H"
 #include "labelHashSet.H"
 #include "Map.H"
+#include "EdgeMap.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -573,6 +574,14 @@ public:
                     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
                 (
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
index 041b831ac19..3106643bb00 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C
@@ -1069,6 +1069,147 @@ bool 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 primitiveMesh::checkEdgeAlignment
+(
+    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;
+    }
+
+    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);
+        }
+    }
+
+    if (nDirs == vector::nComponents)
+    {
+        return false;
+    }
+
+
+
+    const pointField& p = points();
+    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 primitiveMesh::checkUpperTriangular
 (
     const bool report,
-- 
GitLab