reconstructParMesh: Match face point averages on coupled patches

In the event that matching centroids across a coupled patch pair fails,
we fall back to matching the face point average. The latter can be
obtained more reliably on degenerate faces as the calculation does not
involve division by the face area.

This fallback was already implemented as part of processorPolyPatch.
This change also applies it to the faceCoupleInfo class used by
reconstructParMesh.

src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1475,12 +1475,13 @@ void Foam::faceCoupleInfo::perfectPointMatch
         // Faces do not have to be ordered (but all have
         // to match). Note: Faces will be already ordered if we enter here from
         // construct from meshes.
+
         matchedAllFaces = matchPoints
         (
             calcFaceCentres<List>
             (
                 cutFaces(),
-                cutPoints_,
+                cutFaces().points(),
                 0,
                 cutFaces().size()
             ),
@@ -1492,9 +1493,43 @@ void Foam::faceCoupleInfo::perfectPointMatch
                 slavePatch().size()
             ),
             scalarField(slavePatch().size(), absTol),
-            true,
+            false,
             cutToSlaveFaces_
         );
+
+        // If some of the face centres did not match, then try to match the
+        // point averages instead. There is no division by the face area in
+        // calculating the point average, so this is more stable when faces
+        // collapse onto a line or point.
+        if (!matchedAllFaces)
+        {
+            labelList cutToSlaveFacesTemp(cutToSlaveFaces_.size(), -1);
+
+            matchPoints
+            (
+                calcFacePointAverages<List>
+                (
+                    cutFaces(),
+                    cutFaces().points(),
+                    0,
+                    cutFaces().size()
+                ),
+                calcFacePointAverages<IndirectList>
+                (
+                    slavePatch(),
+                    slavePatch().points(),
+                    0,
+                    slavePatch().size()
+                ),
+                scalarField(slavePatch().size(), absTol),
+                true,
+                cutToSlaveFacesTemp
+            );
+
+            cutToSlaveFaces_ = max(cutToSlaveFaces_, cutToSlaveFacesTemp);
+
+            matchedAllFaces = min(cutToSlaveFaces_) != -1;
+        }
     }
 
 

src/dynamicMesh/polyMeshAdder/faceCoupleInfo.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -218,6 +218,16 @@ class faceCoupleInfo
                 const label size
             );
 
+            //- Calculate face point averages from (subset of) faces.
+            template<template<class> class FaceList>
+            static pointField calcFacePointAverages
+            (
+                const FaceList<face>&,
+                const pointField&,
+                const label start,
+                const label size
+            );
+
             //- Write edges
             static void writeOBJ
             (

src/dynamicMesh/polyMeshAdder/faceCoupleInfoTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -49,4 +49,29 @@ Foam::pointField Foam::faceCoupleInfo::calcFaceCentres
 }
 
 
+template<template<class> class FaceList>
+Foam::pointField Foam::faceCoupleInfo::calcFacePointAverages
+(
+    const FaceList<face>& faces,
+    const pointField& points,
+    const label start,
+    const label size
+)
+{
+    pointField fpa(size, Zero);
+
+    label facei = start;
+
+    forAll(fpa, i)
+    {
+        forAll(faces[facei], j)
+        {
+            fpa[i] += points[faces[facei][j]];
+        }
+        fpa[i] /= faces[facei++].size();
+    }
+    return fpa;
+}
+
+
 // ************************************************************************* //