From f7cbb0672beada63b88806ed62e4665868a23465 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Thu, 10 Mar 2022 14:32:19 +0100
Subject: [PATCH] ENH: globalIndex gather for PatchTools::gatherAndMerge
 (#2402)

- lower memory overhead, simpler code and eliminates need for
  ListListOps::combineOffset()

- optional handling of local faces/points for re-using in different
  contexts

STYLE: labelUList instead of labelList for globalMesh mergePoints

STYLE: adjust verbose information from mergePoints

- also report the current new-point location
---
 src/OpenFOAM/meshes/meshTools/mergePoints.C   |  47 ++++---
 .../polyMesh/globalMeshData/globalMeshData.C  |   8 +-
 .../polyMesh/globalMeshData/globalMeshData.H  |   2 +-
 .../primitiveMesh/PatchTools/PatchTools.H     |  15 ++-
 .../PatchTools/PatchToolsGatherAndMerge.C     |  83 ++++++------
 .../PrimitivePatch/PrimitivePatch.H           |   6 +
 .../surfaceFieldValue/surfaceFieldValue.C     | 119 +++++-------------
 7 files changed, 129 insertions(+), 151 deletions(-)

diff --git a/src/OpenFOAM/meshes/meshTools/mergePoints.C b/src/OpenFOAM/meshes/meshTools/mergePoints.C
index e095952f24..b5d6049755 100644
--- a/src/OpenFOAM/meshes/meshTools/mergePoints.C
+++ b/src/OpenFOAM/meshes/meshTools/mergePoints.C
@@ -89,12 +89,12 @@ Foam::label Foam::mergePoints
 
 
     Field<scalar> sortedTol(nPoints);
-    forAll(order, sortI)
+    forAll(order, sorti)
     {
-        const point_type& pt = points[order[sortI]];
+        const point_type& pt = points[order[sorti]];
 
         // Use scalar precision
-        sortedTol[sortI] =
+        sortedTol[sorti] =
             2*mergeTol*
             (
                 mag(scalar(pt.x() - compareOrigin.x()))
@@ -109,11 +109,16 @@ Foam::label Foam::mergePoints
     label pointi = order[0];
     pointMap[pointi] = newPointi++;
 
-    for (label sortI = 1; sortI < order.size(); ++sortI)
+    /// if (verbose)
+    /// {
+    ///     Pout<< "Foam::mergePoints : [0] Uniq point " << pointi << endl;
+    /// }
+
+    for (label sorti = 1; sorti < order.size(); ++sorti)
     {
         // Get original point index
-        const label pointi = order[sortI];
-        const scalar mag2 = magSqrDist[order[sortI]];
+        const label pointi = order[sorti];
+        const scalar mag2 = magSqrDist[order[sorti]];
 
         // Convert to scalar precision
         // NOTE: not yet using point_type template parameter
@@ -130,13 +135,13 @@ Foam::label Foam::mergePoints
 
         for
         (
-            label prevSortI = sortI - 1;
-            prevSortI >= 0
-         && (mag(magSqrDist[order[prevSortI]] - mag2) <= sortedTol[sortI]);
-            --prevSortI
+            label prevSorti = sorti - 1;
+            prevSorti >= 0
+         && (mag(magSqrDist[order[prevSorti]] - mag2) <= sortedTol[sorti]);
+            --prevSorti
         )
         {
-            const label prevPointi = order[prevSortI];
+            const label prevPointi = order[prevSorti];
 
             // Convert to scalar precision
             // NOTE: not yet using point_type template parameter
@@ -164,16 +169,22 @@ Foam::label Foam::mergePoints
 
             if (verbose)
             {
-                Pout<< "Foam::mergePoints : Merging points "
-                    << pointi << " and " << equalPointi
-                    << " with coordinates:" << points[pointi]
-                    << " and " << points[equalPointi]
-                    << endl;
+                Pout<< "Foam::mergePoints : [" << pointMap[pointi]
+                    << "] Point " << pointi << " duplicate of " << equalPointi
+                    << " : coordinates:" << points[pointi]
+                    << " and " << points[equalPointi] << endl;
             }
         }
         else
         {
             // Differs. Store new point.
+
+            /// if (verbose)
+            /// {
+            ///     Pout<< "Foam::mergePoints : [" << newPointi
+            ///         << "] Uniq point " << pointi << endl;
+            /// }
+
             pointMap[pointi] = newPointi++;
         }
     }
@@ -200,7 +211,7 @@ bool Foam::mergePoints
     typename PointList::const_reference origin
 )
 {
-    const label nUnique = mergePoints
+    const label nUnique = Foam::mergePoints
     (
         points,
         mergeTol,
@@ -209,7 +220,7 @@ bool Foam::mergePoints
         origin
     );
 
-    newPoints.setSize(nUnique);
+    newPoints.resize_nocopy(nUnique);
     forAll(pointMap, pointi)
     {
         newPoints[pointMap[pointi]] = points[pointi];
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index d08c7be033..a16d392b28 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -138,7 +138,7 @@ void Foam::globalMeshData::calcSharedPoints() const
 
     // Calculate all shared points (exclude points that are only
     // on two coupled patches). This does all the hard work.
-    globalPoints parallelPoints(mesh_, false, true);
+    const globalPoints parallelPoints(mesh_, false, true);
 
     // Count the number of master points
     label nMaster = 0;
@@ -155,7 +155,7 @@ void Foam::globalMeshData::calcSharedPoints() const
     }
 
     // Allocate global numbers
-    globalIndex masterNumbering(nMaster);
+    const globalIndex masterNumbering(nMaster);
 
     nGlobalPoints_ = masterNumbering.totalSize();
 
@@ -2472,7 +2472,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
 
 Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
 (
-    const labelList& meshPoints,
+    const labelUList& meshPoints,
     const Map<label>& /* unused: meshPointMap */,
     labelList& pointToGlobal,
     labelList& uniqueMeshPoints
@@ -2498,7 +2498,7 @@ Foam::autoPtr<Foam::globalIndex> Foam::globalMeshData::mergePoints
     // - from coupled point to global patch point
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    globalIndex globalPPoints(meshPoints.size());
+    const globalIndex globalPPoints(meshPoints.size());
 
     labelList patchToCoupled(meshPoints.size(), -1);
     label nCoupled = 0;
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
index 8a58b914aa..41ed3e36fb 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.H
@@ -576,7 +576,7 @@ public:
                 //  - the global number for all local points.
                 autoPtr<globalIndex> mergePoints
                 (
-                    const labelList& meshPoints,
+                    const labelUList& meshPoints,
                     const Map<label>& meshPointMap,  //!< currently unused
                     labelList& pointToGlobal,
                     labelList& uniqueMeshPoints
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
index a50b0fdf32..5272bf8069 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
@@ -218,12 +218,20 @@ public:
 
 
     //- Gather points and faces onto master and merge into single patch.
-    //  Note: uses faces/points, not localFaces/localPoints.
+    //  Note: Normally uses faces/points (not localFaces/localPoints)
+    //
+    //  \param[in] mergeDist Geometric merge tolerance for Foam::mergePoints
+    //  \param[in] pp The patch to merge
+    //  \param[out] mergedPoints
+    //  \param[out] mergedFaces
+    //  \param[out] pointMergeMap
+    //  \param[in] useLocal gather/merge the localFaces/localPoints
+    //      instead of (global) faces/points
     template<class FaceList, class PointField>
     static void gatherAndMerge
     (
         const scalar mergeDist,
-        const PrimitivePatch<FaceList, PointField>& p,
+        const PrimitivePatch<FaceList, PointField>& pp,
         Field
         <
             typename PrimitivePatch<FaceList, PointField>::point_type
@@ -232,7 +240,8 @@ public:
         <
             typename PrimitivePatch<FaceList, PointField>::face_type
         >& mergedFaces,
-        labelList& pointMergeMap
+        labelList& pointMergeMap,
+        const bool useLocal = false
     );
 
     //- Gather (mesh!) points and faces onto master and merge collocated
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
index 6e22d5140a..305673e989 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
@@ -37,7 +37,7 @@ template<class FaceList, class PointField>
 void Foam::PatchTools::gatherAndMerge
 (
     const scalar mergeDist,
-    const PrimitivePatch<FaceList, PointField>& p,
+    const PrimitivePatch<FaceList, PointField>& pp,
     Field
     <
         typename PrimitivePatch<FaceList, PointField>::point_type
@@ -46,73 +46,86 @@ void Foam::PatchTools::gatherAndMerge
     <
         typename PrimitivePatch<FaceList, PointField>::face_type
     >& mergedFaces,
-    labelList& pointMergeMap
+    labelList& pointMergeMap,
+    const bool useLocal
 )
 {
     typedef typename PrimitivePatch<FaceList,PointField>::face_type FaceType;
     typedef typename PrimitivePatch<FaceList,PointField>::point_type PointType;
 
-    // Collect points from all processors
-    labelList pointSizes;
+    // Faces from all ranks
+    const globalIndex faceAddr(pp.size(), globalIndex::gatherOnly{});
+
+    // Points from all ranks
+    const globalIndex pointAddr
+    (
+        (useLocal ? pp.localPoints().size() : pp.points().size()),
+        globalIndex::gatherOnly{}
+    );
+
+    if (useLocal)
     {
-        const globalIndex gi(p.points().size());
-
-        gi.gather(p.points(), mergedPoints);
-
-        pointSizes = gi.sizes();
+        faceAddr.gather(pp.localFaces(), mergedFaces);
+        pointAddr.gather(pp.localPoints(), mergedPoints);
+    }
+    else
+    {
+        faceAddr.gather(pp, mergedFaces);
+        pointAddr.gather(pp.points(), mergedPoints);
     }
 
-    // Collect faces from all processors and renumber using sizes of
-    // gathered points
+    // Relabel faces according to global point offsets
+    for (const label proci : faceAddr.subProcs())
     {
-        List<List<FaceType>> gatheredFaces(Pstream::nProcs());
-        gatheredFaces[Pstream::myProcNo()] = p;
-        Pstream::gatherList(gatheredFaces);
+        SubList<FaceType> slot(mergedFaces, faceAddr.range(proci));
 
-        if (Pstream::master())
+        for (auto& f : slot)
         {
-            mergedFaces = static_cast<const List<FaceType>&>
-            (
-                ListListOps::combineOffset<List<FaceType>>
-                (
-                    gatheredFaces,
-                    pointSizes,
-                    accessOp<List<FaceType>>(),
-                    offsetOp<FaceType>()
-                )
-            );
+            pointAddr.inplaceToGlobal(proci, f);
         }
     }
 
-    if (Pstream::master())
+    // Merging points
+    bool hasMerged = false;
+
+    if (Pstream::parRun() && Pstream::master())
     {
         Field<PointType> newPoints;
         labelList oldToNew;
 
-        bool hasMerged = mergePoints
+        hasMerged = Foam::mergePoints
         (
             mergedPoints,
             mergeDist,
-            false,                  // verbosity
+            false,  // verbosity
             oldToNew,
             newPoints
         );
 
         if (hasMerged)
         {
-            // Store point mapping
-            pointMergeMap.transfer(oldToNew);
-
-            // Copy points
-            mergedPoints.transfer(newPoints);
-
             // Relabel faces
             for (auto& f : mergedFaces)
             {
-                inplaceRenumber(pointMergeMap, f);
+                inplaceRenumber(oldToNew, f);
+            }
+
+            // Store newly merged points
+            mergedPoints.transfer(newPoints);
+
+            // Store point mapping
+            if (notNull(pointMergeMap))
+            {
+                pointMergeMap.transfer(oldToNew);
             }
         }
     }
+
+    // TDB:
+    // if (!hasMerged && notNull(pointMergeMap))
+    // {
+    //     pointMergeMap.clear();
+    // }
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
index 13e198de3a..a9f47301d1 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
@@ -468,9 +468,15 @@ public:
             bool hasPointNormals() const { return bool(pointNormalsPtr_); }
 
             bool hasBoundaryPoints() const { return bool(boundaryPointsPtr_); }
+
+            // These ones are currently all calculated together:
+            // - edges(), faceFaces(), edgeFaces(), faceEdges()
+
+            bool hasEdges() const { return bool(edgesPtr_); }
             bool hasFaceFaces() const { return bool(faceFacesPtr_); }
             bool hasEdgeFaces() const { return bool(edgeFacesPtr_); }
             bool hasFaceEdges() const { return bool(faceEdgesPtr_); }
+
             bool hasPointEdges() const { return bool(pointEdgesPtr_); }
             bool hasPointFaces() const { return bool(pointFacesPtr_); }
 
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
index 92a5839027..74480573fe 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +32,7 @@ License
 #include "coupledPolyPatch.H"
 #include "sampledSurface.H"
 #include "mergePoints.H"
-#include "uindirectPrimitivePatch.H"
+#include "indirectPrimitivePatch.H"
 #include "PatchTools.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -353,104 +353,43 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
     pointField& points
 ) const
 {
-    List<faceList> allFaces(Pstream::nProcs());
-    List<pointField> allPoints(Pstream::nProcs());
+    labelList whichFaces(faceId_);
 
+    // Remap patch-face ids to mesh face ids
+    forAll(whichFaces, i)
     {
-        IndirectList<face> selectedFaces(mesh_.faces(), labelList(faceId_));
-        labelList& meshFaceIds = selectedFaces.addressing();
-
-        forAll(meshFaceIds, i)
+        const label patchi = facePatchId_[i];
+        if (patchi != -1)
         {
-            const label patchi = facePatchId_[i];
-            if (patchi != -1)
-            {
-                meshFaceIds[i] += mesh_.boundaryMesh()[patchi].start();
-            }
-        }
-
-        // Add local faces and points to the all* lists
-        uindirectPrimitivePatch pp(selectedFaces, mesh_.points());
-
-        allFaces[Pstream::myProcNo()] = pp.localFaces();
-        allPoints[Pstream::myProcNo()] = pp.localPoints();
-    }
-
-    Pstream::gatherList(allFaces);
-    Pstream::gatherList(allPoints);
-
-    // Renumber and flatten
-    label nFaces = 0;
-    label nPoints = 0;
-    forAll(allFaces, proci)
-    {
-        nFaces += allFaces[proci].size();
-        nPoints += allPoints[proci].size();
-    }
-
-    faces.resize(nFaces);
-    points.resize(nPoints);
-
-    nFaces = 0;
-    nPoints = 0;
-
-    // My data first
-    {
-        for (const face& f : allFaces[Pstream::myProcNo()])
-        {
-            faces[nFaces++] = offsetOp<face>()(f, nPoints);
-        }
-
-        for (const point& p : allPoints[Pstream::myProcNo()])
-        {
-            points[nPoints++] = p;
+            whichFaces[i] += mesh_.boundaryMesh()[patchi].start();
         }
     }
 
-    // Other proc data follows
-    forAll(allFaces, proci)
-    {
-        if (proci == Pstream::myProcNo())
-        {
-            continue;
-        }
-
-        for (const face& f : allFaces[proci])
-        {
-            faces[nFaces++] = offsetOp<face>()(f, nPoints);
-        }
-
-        for (const point& p : allPoints[proci])
-        {
-            points[nPoints++] = p;
-        }
-    }
-
-    // Merge
-    labelList oldToNew;
-    pointField newPoints;
-    bool hasMerged = mergePoints
+    indirectPrimitivePatch pp
     (
-        points,
-        SMALL,
-        false,
-        oldToNew,
-        newPoints
+        IndirectList<face>(mesh_.faces(), std::move(whichFaces)),
+        mesh_.points()
     );
 
-    if (hasMerged)
-    {
-        if (debug)
-        {
-            Pout<< "Merged from " << points.size()
-                << " down to " << newPoints.size() << " points" << endl;
-        }
 
-        points.transfer(newPoints);
-        for (face& f : faces)
-        {
-            inplaceRenumber(oldToNew, f);
-        }
+    if (Pstream::parRun())
+    {
+        labelList pointsMap;
+
+        PatchTools::gatherAndMerge
+        (
+            SMALL,  // mergeDist
+            pp,
+            points,
+            faces,
+            pointsMap,
+            true   // useLocal=true
+        );
+    }
+    else
+    {
+        faces = pp.localFaces();
+        points = pp.localPoints();
     }
 }