From a242a3787b474dccbe8f8c1b030e189cd56e6285 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 24 Aug 2011 08:58:48 +0100
Subject: [PATCH] ENH: checkMesh: added check for coupled point synchronisation

---
 .../manipulation/checkMesh/checkGeometry.C    | 139 +++++++++++++++++-
 .../manipulation/checkMesh/checkGeometry.H    |   4 +
 .../basic/coupled/coupledPolyPatch.H          |  19 +--
 3 files changed, 152 insertions(+), 10 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
index 82e0dedd82..9004ddade9 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C
@@ -262,6 +262,125 @@ bool Foam::checkWedges
 }
 
 
+bool Foam::checkCoupledPoints
+(
+    const polyMesh& mesh,
+    const bool report,
+    labelHashSet* setPtr
+)
+{
+    const pointField& p = mesh.points();
+    const faceList& fcs = mesh.faces();
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Zero'th point on coupled faces
+    pointField nbrZeroPoint(fcs.size()-mesh.nInternalFaces(), vector::max);
+
+    // Exchange zero point
+    forAll(patches, patchI)
+    {
+        if (patches[patchI].coupled())
+        {
+            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+            (
+                patches[patchI]
+            );
+
+            forAll(cpp, i)
+            {
+                label bFaceI = cpp.start()+i-mesh.nInternalFaces();
+                const point& p0 = p[cpp[i][0]];
+                nbrZeroPoint[bFaceI] = p0;
+            }
+        }
+    }
+    syncTools::swapBoundaryFacePositions(mesh, nbrZeroPoint);
+
+    // Compare to local ones. Use same tolerance as for matching
+    label nErrorFaces = 0;
+    scalar avgMismatch = 0;
+    label nCoupledFaces = 0;
+
+    forAll(patches, patchI)
+    {
+        if (patches[patchI].coupled())
+        {
+            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+            (
+                patches[patchI]
+            );
+
+            if (cpp.owner())
+            {
+                scalarField smallDist
+                (
+                    cpp.calcFaceTol
+                    (
+                        cpp.matchTolerance(),
+                        cpp,
+                        cpp.points(),
+                        cpp.faceCentres()
+                    )
+                );
+
+                forAll(cpp, i)
+                {
+                    label bFaceI = cpp.start()+i-mesh.nInternalFaces();
+                    const point& p0 = p[cpp[i][0]];
+
+                    scalar d = mag(p0 - nbrZeroPoint[bFaceI]);
+
+                    if (d > smallDist[i])
+                    {
+                        if (setPtr)
+                        {
+                            setPtr->insert(cpp.start()+i);
+                        }
+                        nErrorFaces++;
+                    }
+                    avgMismatch += d;
+                    nCoupledFaces++;
+                }
+            }
+        }
+    }
+
+    reduce(nErrorFaces, sumOp<label>());
+    reduce(avgMismatch, maxOp<scalar>());
+    reduce(nCoupledFaces, sumOp<label>());
+
+    if (nCoupledFaces > 0)
+    {
+        avgMismatch /= nCoupledFaces;
+    }
+
+    if (nErrorFaces > 0)
+    {
+        if (report)
+        {
+            Info<< "  **Error in coupled point location: "
+                << nErrorFaces
+                << " faces have their 0th vertex not opposite"
+                << " their coupled equivalent. Average mismatch "
+                << avgMismatch << "."
+                << endl;
+        }
+
+        return true;
+    }
+    else
+    {
+        if (report)
+        {
+            Info<< "    Coupled point location match (average "
+                << avgMismatch << ") OK." << endl;
+        }
+
+        return false;
+    }
+}
+
+
 Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
 {
     label noFailedChecks = 0;
@@ -451,6 +570,25 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
+    {
+        faceSet faces(mesh, "coupledFaces", mesh.nFaces()/100 + 1);
+        if (checkCoupledPoints(mesh, true, &faces))
+        {
+            noFailedChecks++;
+
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            if (nFaces > 0)
+            {
+                Info<< "  <<Writing " << nFaces
+                    << " faces with incorrectly matched 0th vertex to set "
+                    << faces.name() << endl;
+                faces.instance() = mesh.pointsInstance();
+                faces.write();
+            }
+        }
+    }
+
     if (allGeometry)
     {
         faceSet faces(mesh, "lowQualityTetFaces", mesh.nFaces()/100+1);
@@ -591,6 +729,5 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
-
     return noFailedChecks;
 }
diff --git a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H
index 8392e551ea..edc1f44826 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H
+++ b/applications/utilities/mesh/manipulation/checkMesh/checkGeometry.H
@@ -9,6 +9,7 @@ namespace Foam
 
     label findOppositeWedge(const polyMesh&, const wedgePolyPatch&);
 
+    //- Check wedge orientation
     bool checkWedges
     (
         const polyMesh&,
@@ -17,5 +18,8 @@ namespace Foam
         labelHashSet*
     );
 
+    //- Check 0th vertex on coupled faces
+    bool checkCoupledPoints(const polyMesh&, const bool report, labelHashSet*);
+
     label checkGeometry(const polyMesh& mesh, const bool allGeometry);
 }
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
index 0514db1ab8..b3368d5cfb 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H
@@ -156,15 +156,6 @@ protected:
             const pointField&
         );
 
-        //- Calculate typical tolerance per face. Is currently max distance
-        //  from face centre to any of the face vertices.
-        static scalarField calcFaceTol
-        (
-            const UList<face>& faces,
-            const pointField& points,
-            const pointField& faceCentres
-        );
-
         //- Get the number of vertices face f needs to be rotated such that
         //  its f[0] point aligns with given anchor (within tol).
         static label getRotation
@@ -332,6 +323,16 @@ public:
             labelList& rotation
         ) const = 0;
 
+        //- Calculate typical tolerance per face. Is currently max distance
+        //  from face centre to any of the face vertices.
+        static scalarField calcFaceTol
+        (
+            const scalar matchTol,
+            const UList<face>& faces,
+            const pointField& points,
+            const pointField& faceCentres
+        );
+
         //- Write the polyPatch data as a dictionary
         virtual void write(Ostream&) const;
 };