From b37e628ac8001d5b3e44390c54e234e34cfafac6 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Mon, 6 Nov 2017 14:34:49 +0000
Subject: [PATCH] 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.
---
 .../polyMeshAdder/faceCoupleInfo.C            | 41 +++++++++++++++++--
 .../polyMeshAdder/faceCoupleInfo.H            | 12 +++++-
 .../polyMeshAdder/faceCoupleInfoTemplates.C   | 27 +++++++++++-
 3 files changed, 75 insertions(+), 5 deletions(-)

diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
index b4f8f7449f..3b691ba0f5 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.C
+++ b/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;
+        }
     }
 
 
diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.H b/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.H
index b5d7b1fe28..28137d41b1 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfo.H
+++ b/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
             (
diff --git a/src/dynamicMesh/polyMeshAdder/faceCoupleInfoTemplates.C b/src/dynamicMesh/polyMeshAdder/faceCoupleInfoTemplates.C
index 3b25a6da51..2158c0c92a 100644
--- a/src/dynamicMesh/polyMeshAdder/faceCoupleInfoTemplates.C
+++ b/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;
+}
+
+
 // ************************************************************************* //