From 24c0b30d48fc4722aa6d40ac4d0d163d15cf14f9 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Tue, 15 Mar 2022 21:38:49 +0100
Subject: [PATCH] ENH: mergePoints and patch gatherAndMerge improvements
 (#2402)

- when writing surface formats (eg, vtk, ensight etc) the sampled
  surfaces merge the faces/points originating from different
  processors into a single surface (ie, patch gatherAndMerge).

  Previous versions of mergePoints simply merged all points possible,
  which proves to be rather slow for larger meshes. This has now been
  modified to only consider boundary points, which reduces the number
  of points to consider. As part of this change, the reference point
  is now always equivalent to the min of the bounding box, which
  reduces the number of search loops. The merged points retain their
  original order.

- inplaceMergePoints version to simplify use and improve code
  robustness and efficiency.

ENH: make PrimitivePatch::boundaryPoints() less costly

- if edge addressing does not already exist, it will now simply walk
  the local face edges directly to define the boundary points.

  This avoids a rather large overhead of the full faceFaces,
  edgeFaces, faceEdges addressing.

  This operation is now more important since it is used in the revised
  patch gatherAndMerge.

ENH: topological merge for mesh-based surfaces in surfaceFieldValue
---
 .../manipulation/polyDualMesh/meshDualiser.C  |   4 +-
 src/OpenFOAM/meshes/meshTools/mergePoints.C   | 515 ++++++++++++++----
 src/OpenFOAM/meshes/meshTools/mergePoints.H   | 217 +++++++-
 .../polyMesh/globalMeshData/globalMeshData.C  |  12 +-
 .../primitiveMesh/PatchTools/PatchTools.H     |  15 +-
 .../PatchTools/PatchToolsGatherAndMerge.C     |  62 ++-
 .../PrimitivePatch/PrimitivePatch.H           |   8 +-
 .../PrimitivePatch/PrimitivePatchAddressing.C |   4 +-
 .../PrimitivePatch/PrimitivePatchBdryPoints.C |  77 ++-
 src/dynamicMesh/polyMeshAdder/polyMeshAdder.C |   6 +-
 .../surfaceFieldValue/surfaceFieldValue.C     |  26 +-
 .../snappyHexMeshDriver/snappyLayerDriver.C   |   4 +-
 .../snappyHexMeshDriver/snappySnapDriver.C    |   2 +-
 .../advancingFrontAMIParallelOps.C            |  13 +-
 src/meshTools/edgeMesh/edgeMesh.C             |  15 +-
 .../extendedEdgeMesh/extendedEdgeMesh.C       |   2 +-
 .../surfaceIntersection/surfaceIntersection.C |  12 +-
 .../triSurfaceTools/triSurfaceTools.C         |  17 +-
 .../distributedTriSurfaceMesh.C               |   6 +-
 .../meshToMesh/meshToMeshParallelOps.C        |  17 +-
 .../sampledSet/patchEdge/patchEdgeSet.C       |   9 +-
 .../surface/isoSurface/isoSurfaceCell.C       |  43 +-
 .../surface/isoSurface/isoSurfacePoint.C      |   6 +-
 src/surfMesh/MeshedSurface/MeshedSurface.C    |  17 +-
 src/surfMesh/surfaceFormats/tri/TRIReader.C   |   1 +
 src/surfMesh/triSurface/triSurfaceStitch.C    |  19 +-
 26 files changed, 833 insertions(+), 296 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C b/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C
index c512857b76..5fcf5d2f89 100644
--- a/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C
+++ b/applications/utilities/mesh/manipulation/polyDualMesh/meshDualiser.C
@@ -55,7 +55,7 @@ void Foam::meshDualiser::checkPolyTopoChange(const polyTopoChange& meshMod)
     }
 
     labelList oldToNew;
-    label nUnique = mergePoints
+    label nUnique = Foam::mergePoints
     (
         points,
         1e-6,
@@ -223,7 +223,7 @@ Foam::label Foam::meshDualiser::addInternalFace
         pointField facePoints(meshMod.points(), newFace);
 
         labelList oldToNew;
-        label nUnique = mergePoints
+        label nUnique = Foam::mergePoints
         (
             facePoints,
             1e-6,
diff --git a/src/OpenFOAM/meshes/meshTools/mergePoints.C b/src/OpenFOAM/meshes/meshTools/mergePoints.C
index b5d6049755..56746484a8 100644
--- a/src/OpenFOAM/meshes/meshTools/mergePoints.C
+++ b/src/OpenFOAM/meshes/meshTools/mergePoints.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2019 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,49 +26,47 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "ListOps.H"
-#include "point.H"
-#include "Field.H"
+#include "ListOps.H"  // sortedOrder, ListOps::identity
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Implementation  * * * * * * * * * * * * * * //
 
-template<class PointList>
-Foam::label Foam::mergePoints
+template<class PointList, class IndexerOp>
+Foam::label Foam::Detail::mergePoints
 (
     const PointList& points,
+    const IndexerOp& indexer,
+    const label nSubPoints,
+    labelList& pointToUnique,
+    labelList& uniquePoints,
     const scalar mergeTol,
-    const bool verbose,
-    labelList& pointMap,
-    typename PointList::const_reference origin
+    const bool verbose
 )
 {
-    typedef typename PointList::value_type point_type;
+    const label nTotPoints = points.size();
 
-    const label nPoints = points.size();
-
-    // Create an old to new point mapping array
-    pointMap.resize_nocopy(nPoints);
-    pointMap = -1;
-
-    if (!nPoints)
+    if (!nTotPoints || !nSubPoints)
     {
-        return 0;
+        // Nothing to do
+        pointToUnique = identity(nTotPoints);
+        uniquePoints = pointToUnique;
+        return 0;  // No points removed
     }
 
-    point_type compareOrigin = origin;
-    if (origin == point_type::max)
+    // Properly size for old to new mapping array
+    pointToUnique.resize_nocopy(nTotPoints);
+
+
+    // Use the boundBox minimum as the reference point. This
+    // stretches distances with fewer collisions than a mid-point
+    // reference would.
+
+    auto comparePoint(points[indexer(0)]);
+    for (label pointi = 1; pointi < nSubPoints; ++pointi)
     {
-        // Use average of input points to define a comparison origin.
-        // Same as sum(points)/nPoints, but handles different list types
-        compareOrigin = points[0];
-        for (label pointi=1; pointi < nPoints; ++pointi)
-        {
-            compareOrigin += points[pointi];
-        }
-        compareOrigin /= nPoints;
+        comparePoint = min(comparePoint, points[indexer(pointi)]);
     }
 
-    // We're comparing distance squared to origin first.
+    // We're comparing distance squared to reference point first.
     // Say if starting from two close points:
     //     x, y, z
     //     x+mergeTol, y+mergeTol, z+mergeTol
@@ -77,156 +75,445 @@ Foam::label Foam::mergePoints
     //     x^2+y^2+z^2 + 2*mergeTol*(x+z+y) + mergeTol^2*...
     // so the difference will be 2*mergeTol*(x+y+z)
 
-    const scalar mergeTolSqr = Foam::sqr(mergeTol);
+    const scalar mergeTolSqr(magSqr(mergeTol));
 
-    // Sort points by magSqr
-    List<scalar> magSqrDist(nPoints);
-    forAll(points, pointi)
+    // Use magSqr distance for the points sort order
+    List<scalar> sqrDist(nSubPoints);
+    for (label pointi = 0; pointi < nSubPoints; ++pointi)
     {
-        magSqrDist[pointi] = magSqr(points[pointi] - compareOrigin);
+        const auto& p = points[indexer(pointi)];
+
+        sqrDist[pointi] =
+        (
+            // Use scalar precision
+            magSqr(scalar(p.x() - comparePoint.x()))
+          + magSqr(scalar(p.y() - comparePoint.y()))
+          + magSqr(scalar(p.z() - comparePoint.z()))
+        );
     }
-    labelList order(Foam::sortedOrder(magSqrDist));
+    labelList order(Foam::sortedOrder(sqrDist));
 
-
-    Field<scalar> sortedTol(nPoints);
+    List<scalar> sortedTol(nSubPoints);
     forAll(order, sorti)
     {
-        const point_type& pt = points[order[sorti]];
+        const auto& p = points[indexer(order[sorti])];
 
-        // Use scalar precision
         sortedTol[sorti] =
+        (
             2*mergeTol*
             (
-                mag(scalar(pt.x() - compareOrigin.x()))
-              + mag(scalar(pt.y() - compareOrigin.y()))
-              + mag(scalar(pt.z() - compareOrigin.z()))
-            );
+                // Use scalar precision
+                mag(scalar(p.x() - comparePoint.x()))
+              + mag(scalar(p.y() - comparePoint.y()))
+              + mag(scalar(p.z() - comparePoint.z()))
+            )
+        );
     }
 
-    label newPointi = 0;
 
-    // Handle 0th point separately (is always unique)
-    label pointi = order[0];
-    pointMap[pointi] = newPointi++;
+    // Bookkeeping parameters
+    // ~~~~~~~~~~~~~~~~~~~~~~
 
-    /// if (verbose)
-    /// {
-    ///     Pout<< "Foam::mergePoints : [0] Uniq point " << pointi << endl;
-    /// }
+    // Will only be working on a subset of the points
+    // Can use a slice of pointToUnique (full length not needed until later).
+    SubList<label> subPointMap(pointToUnique, nSubPoints);
 
-    for (label sorti = 1; sorti < order.size(); ++sorti)
+    // Track number of unique points - this will form an offsets table
+    labelList newPointCounts(nSubPoints, Zero);
+
+    label nNewPoints = 0;
+    for (label sorti = 0; sorti < order.size(); ++sorti)
     {
-        // Get original point index
+        // The (sub)point index
         const label pointi = order[sorti];
-        const scalar mag2 = magSqrDist[order[sorti]];
+        const scalar currDist = sqrDist[order[sorti]];
+        const auto& currPoint = points[indexer(pointi)];
 
-        // Convert to scalar precision
-        // NOTE: not yet using point_type template parameter
-        const point pt
-        (
-            scalar(points[pointi].x()),
-            scalar(points[pointi].y()),
-            scalar(points[pointi].z())
-        );
+        // Compare to previous points to find equal one
+        // - automatically a no-op for sorti == 0 (the first point)
 
-
-        // Compare to previous points to find equal one.
-        label equalPointi = -1;
+        bool matched = false;
 
         for
         (
             label prevSorti = sorti - 1;
-            prevSorti >= 0
-         && (mag(magSqrDist[order[prevSorti]] - mag2) <= sortedTol[sorti]);
+            (
+                prevSorti >= 0
+             && (mag(sqrDist[order[prevSorti]] - currDist) <= sortedTol[sorti])
+            );
             --prevSorti
         )
         {
             const label prevPointi = order[prevSorti];
+            const auto& prevPoint = points[indexer(prevPointi)];
 
-            // Convert to scalar precision
-            // NOTE: not yet using point_type template parameter
-            const point prevPt
+            // Matched within tolerance?
+            matched =
             (
-                scalar(points[prevPointi].x()),
-                scalar(points[prevPointi].y()),
-                scalar(points[prevPointi].z())
+                (
+                    // Use scalar precision
+                    magSqr(scalar(currPoint.x() - prevPoint.x()))
+                  + magSqr(scalar(currPoint.y() - prevPoint.y()))
+                  + magSqr(scalar(currPoint.z() - prevPoint.z()))
+                ) <= mergeTolSqr
             );
 
-            if (magSqr(pt - prevPt) <= mergeTolSqr)
+            if (matched)
             {
-                // Found match.
-                equalPointi = prevPointi;
+                // Both pointi and prevPointi have similar coordinates.
+                // Map to the same new point.
+                subPointMap[pointi] = subPointMap[prevPointi];
 
+                if (verbose)
+                {
+                    Pout<< "Foam::mergePoints : [" << subPointMap[pointi]
+                        << "] Point " << pointi << " duplicate of "
+                        << prevPointi << " : coordinates:" << currPoint
+                        << " and " << prevPoint << endl;
+                }
                 break;
             }
         }
 
-
-        if (equalPointi != -1)
-        {
-            // Same coordinate as equalPointi. Map to same new point.
-            pointMap[pointi] = pointMap[equalPointi];
-
-            if (verbose)
-            {
-                Pout<< "Foam::mergePoints : [" << pointMap[pointi]
-                    << "] Point " << pointi << " duplicate of " << equalPointi
-                    << " : coordinates:" << points[pointi]
-                    << " and " << points[equalPointi] << endl;
-            }
-        }
-        else
+        if (!matched)
         {
             // Differs. Store new point.
+            subPointMap[pointi] = nNewPoints++;
 
-            /// if (verbose)
-            /// {
-            ///     Pout<< "Foam::mergePoints : [" << newPointi
-            ///         << "] Uniq point " << pointi << endl;
-            /// }
-
-            pointMap[pointi] = newPointi++;
+            /// Too verbose
+            ///if (verbose)
+            ///{
+            ///    Pout<< "Foam::mergePoints : [" << subPointMap[pointi]
+            ///        << "] Point " << pointi << endl;
+            ///}
         }
+        ++newPointCounts[subPointMap[pointi]];
     }
 
+    const label nDupPoints(nSubPoints - nNewPoints);
+    const label nUniqPoints(nTotPoints - nDupPoints);
+
     if (verbose)
     {
         Pout<< "Foam::mergePoints : "
-            << newPointi << " of " << points.size() << " unique points"
-            << endl;
+            << "Merging removed " << nDupPoints << '/'
+            << nTotPoints << " points" << endl;
     }
 
-    return newPointi;
+    if (!nDupPoints)
+    {
+        // Nothing to do
+        pointToUnique = identity(nTotPoints);
+        uniquePoints = pointToUnique;
+        return 0;  // No points removed
+    }
+
+
+    // The subPointMap now contains a mapping of the sub-selection
+    // to the list of (sorted) merged points.
+    // Get its sort order to bundle according to the merged point target.
+    // This is in effect an adjacent list of graph edges to mapping back
+    // to the merged points, but in compact form.
+    // Use the previously obtained newPointCounts for the offsets list.
+
+    labelList lookupMerged(std::move(order));
+    Foam::sortedOrder(subPointMap, lookupMerged);
+
+    // Remap inplace to the original points ids
+    for (label& idx : lookupMerged)
+    {
+        idx = indexer(idx);
+    }
+    // The subPointMap slice is not needed beyond here
+
+
+    // Setup initial identity +1 mapping for pointToUnique
+    // The +1 allows negatives to mark duplicates
+
+    ListOps::identity(pointToUnique, 1);
+
+    // The newPointCounts is an offsets table that we use to walk
+    // across the adjacency list (lookupMerged), picking the original
+    // point with the lowest id as the one to retain (master).
+    {
+        label beg = 0;
+        for (const label len : newPointCounts)
+        {
+            if (!len) continue;  // Can be empty
+
+            const label end = (beg + len);
+
+            // Pass 1:
+            // Find the 'master' (lowest point id)
+
+            label masterPointi = lookupMerged[beg];
+
+            for (label iter = beg + 1; iter < end; ++iter)
+            {
+                const label origPointi = lookupMerged[iter];
+
+                if (masterPointi > origPointi)
+                {
+                    masterPointi = origPointi;
+                }
+            }
+
+            // Pass 2:
+            // Markup duplicate points, encoding information about master
+            for (label iter = beg; iter < end; ++iter)
+            {
+                const label origPointi = lookupMerged[iter];
+
+                if (masterPointi != origPointi)
+                {
+                    // Encode the originating 'master' point
+                    pointToUnique[origPointi] = (-masterPointi-1);
+                }
+            }
+
+            beg = end;
+        }
+    }
+
+    // Now have all the information needed
+
+    uniquePoints.resize_nocopy(nUniqPoints);
+    {
+        label uniquei = 0;
+
+        forAll(pointToUnique, pointi)
+        {
+            const label origPointi = pointToUnique[pointi];
+
+            if (origPointi > 0)
+            {
+                // Subtract one to align addressing
+                uniquePoints[uniquei] = (origPointi - 1);
+                pointToUnique[pointi] = uniquei;
+                ++uniquei;
+            }
+            else
+            {
+                // A duplicate point. Also guaranteed that the 'master' point
+                // has a lower index and thus already been seen.
+                const label masterPointi = mag(origPointi) - 1;
+                pointToUnique[pointi] = pointToUnique[masterPointi];
+            }
+        }
+    }
+
+    return nDupPoints;
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class PointList>
+Foam::label Foam::mergePoints
+(
+    const PointList& points,
+    labelList& pointToUnique,
+    labelList& uniquePoints,
+    const scalar mergeTol,
+    const bool verbose
+)
+{
+    const label nTotPoints = points.size();
+
+    if (!nTotPoints)
+    {
+        // Nothing to do
+        pointToUnique.clear();
+        uniquePoints.clear();
+        return 0;  // No points removed
+    }
+
+    return Foam::Detail::mergePoints
+    (
+        points,
+        identityOp(),   // identity indexer
+        nTotPoints,     // == nSubPoints
+        pointToUnique,
+        uniquePoints,
+        mergeTol,
+        verbose
+    );
 }
 
 
 template<class PointList>
-bool Foam::mergePoints
+Foam::label Foam::mergePoints
+(
+    const PointList& points,
+    const labelUList& selection,
+    labelList& pointToUnique,
+    labelList& uniquePoints,
+    const scalar mergeTol,
+    const bool verbose
+)
+{
+    const label nTotPoints = points.size();
+    const label nSubPoints = selection.size();
+
+    if (!nTotPoints || !nSubPoints)
+    {
+        // Nothing to do
+        pointToUnique.clear();
+        uniquePoints.clear();
+        return 0;  // No points removed
+    }
+
+    const auto indexer = [&](const label i) -> label { return selection[i]; };
+
+    return Foam::Detail::mergePoints
+    (
+        points,
+        indexer,
+        nSubPoints,
+        pointToUnique,
+        uniquePoints,
+        mergeTol,
+        verbose
+    );
+}
+
+
+template<class PointList>
+Foam::label Foam::mergePoints
 (
     const PointList& points,
     const scalar mergeTol,
     const bool verbose,
-    labelList& pointMap,
-    List<typename PointList::value_type>& newPoints,
-    typename PointList::const_reference origin
+    labelList& pointToUnique
 )
 {
-    const label nUnique = Foam::mergePoints
+    labelList uniquePoints;
+    const label nChanged = Foam::mergePoints
     (
         points,
+        pointToUnique,
+        uniquePoints,
         mergeTol,
-        verbose,
-        pointMap,
-        origin
+        verbose
     );
 
-    newPoints.resize_nocopy(nUnique);
-    forAll(pointMap, pointi)
+    // Number of unique points
+    return (points.size() - nChanged);
+}
+
+
+template<class PointList>
+Foam::label Foam::inplaceMergePoints
+(
+    PointList& points,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointToUnique
+)
+{
+    labelList uniquePoints;
+    const label nChanged = Foam::mergePoints
+    (
+        points,
+        pointToUnique,
+        uniquePoints,
+        mergeTol,
+        verbose
+    );
+
+    if (nChanged)
     {
-        newPoints[pointMap[pointi]] = points[pointi];
+        // Overwrite
+        points = List<typename PointList::value_type>(points, uniquePoints);
+    }
+    else
+    {
+        // TDB:
+        // pointToUnique.clear();
     }
 
-    return (nUnique != points.size());
+    return nChanged;
+}
+
+
+template<class PointList>
+Foam::label Foam::inplaceMergePoints
+(
+    PointList& points,
+    const labelUList& selection,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointToUnique
+)
+{
+    labelList uniquePoints;
+    const label nChanged = Foam::mergePoints
+    (
+        points,
+        selection,
+        pointToUnique,
+        uniquePoints,
+        mergeTol,
+        verbose
+    );
+
+    if (nChanged)
+    {
+        // Overwrite
+        points = List<typename PointList::value_type>(points, uniquePoints);
+    }
+    else
+    {
+        // TDB:
+        // pointToUnique.clear();
+    }
+
+    return nChanged;
+}
+
+
+template<class PointList>
+Foam::label Foam::mergePoints
+(
+    const PointList& points,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointToUnique,
+    List<typename PointList::value_type>& newPoints
+)
+{
+    const label nTotPoints = points.size();
+
+    if (!nTotPoints)
+    {
+        // Nothing to do
+        pointToUnique.clear();
+        newPoints.clear();
+        return 0;  // No points removed
+    }
+
+    labelList uniquePoints;
+    const label nChanged = Foam::mergePoints
+    (
+        points,
+        pointToUnique,
+        uniquePoints,
+        mergeTol,
+        verbose
+    );
+
+    if (nChanged)
+    {
+        newPoints = List<typename PointList::value_type>(points, uniquePoints);
+    }
+    else
+    {
+        // TDB:
+        // pointToUnique.clear();
+        newPoints = points;
+    }
+
+    return nChanged;
 }
 
 
diff --git a/src/OpenFOAM/meshes/meshTools/mergePoints.H b/src/OpenFOAM/meshes/meshTools/mergePoints.H
index 0954e32412..48330a5f88 100644
--- a/src/OpenFOAM/meshes/meshTools/mergePoints.H
+++ b/src/OpenFOAM/meshes/meshTools/mergePoints.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,15 +25,15 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Description
-    Merge points. See below.
+    Geometric merging of points. See below.
 
 SourceFiles
     mergePoints.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef mergePoints_H
-#define mergePoints_H
+#ifndef Foam_mergePoints_H
+#define Foam_mergePoints_H
 
 #include "scalar.H"
 #include "labelList.H"
@@ -43,35 +43,224 @@ SourceFiles
 namespace Foam
 {
 
+namespace Detail
+{
+
 /*---------------------------------------------------------------------------*\
-                           Function mergePoints Declaration
+                   Function Detail::mergePoints Declaration
 \*---------------------------------------------------------------------------*/
 
-//- Sorts and merges points. All points closer than/equal mergeTol get merged.
-//  Returns the number of unique points and a map from old to new.
+//- Implementation detail for Foam::mergePoints
+//
+//  The implementation algorithm is provided with an \em indexer functor
+//  to map the subset of points to be merged into the full list.
+//  if passed the identityOp this is the same as merging all points.
+//
+//  \returns Number of point duplicates \em removed by merging.
+template<class PointList, class IndexerOp>
+label mergePoints
+(
+    const PointList& points,     //!< The input of all (unmerged) points
+    const IndexerOp& indexer,    //!< Indexer for sub-points into all points
+    const label nSubPoints,      //!< The number of sub-points
+    labelList& pointToUnique,    //!< An old-to-new mapping
+    labelList& uniquePoints,     //!< List of unique points from full list
+    const scalar mergeTol,       //!< Geometric merge tolerance
+    const bool verbose           //!< Debug verbosity
+);
+
+} // End namespace Detail
+
+
+/*---------------------------------------------------------------------------*\
+                       Function mergePoints Declarations
+\*---------------------------------------------------------------------------*/
+
+//- Calculate merge mapping, preserving the original point order.
+//- All points closer/equal mergeTol are to be merged.
+//
+//  \param[in] points The input (unmerged) points
+//  \param[out] pointToUnique  An old-to-new mapping from the original
+//      unmerged point index to the index into unique points.
+//  \param[out] uniquePoints List of unique points from the original
+//      list of unmerged points
+//  \param[in] mergeTol Geometric merge tolerance
+//  \param[in] verbose Debug verbosity
+//
+//  \returns Number of point duplicates \em removed by merging.
+template<class PointList>
+label mergePoints
+(
+    const PointList& points,
+    labelList& pointToUnique,
+    labelList& uniquePoints,
+    const scalar mergeTol = SMALL,
+    const bool verbose = false
+);
+
+
+//- Calculate merge mapping using a subset of points,
+//- preserving the original point order.
+//- All points closer/equal mergeTol are to be merged.
+//
+//  \param[in] allPoints The input of all (unmerged) points
+//  \param[in] selection The subset of points to consider for merging
+//  \param[out] pointToUnique  An old-to-new mapping from the original
+//      unmerged point index to the index into unique points.
+//  \param[out] uniquePoints List of unique points from the original
+//      list of unmerged points
+//  \param[in] mergeTol Geometric merge tolerance
+//  \param[in] verbose Additional debug verbosity
+//
+//  \returns Number of point duplicates \em removed by merging.
+template<class PointList>
+label mergePoints
+(
+    const PointList& points,
+    const labelUList& selection,
+    labelList& pointToUnique,
+    labelList& uniquePoints,
+    const scalar mergeTol = SMALL,
+    const bool verbose = false
+);
+
+
+//- Calculate merge mapping, preserving the original point order.
+//- All points closer/equal mergeTol are to be merged.
+//
+//  \param[in] points The input (unmerged) points
+//  \param[in] mergeTol Geometric merge tolerance
+//  \param[in] verbose Debug verbosity
+//  \param[out] pointToUnique An old-to-new mapping from original
+//      unmerged point index to the index into unique points.
+//
+//  \returns Number of \em unique points after merging.
+//      Can check against the size of points or pointToUnique to determine
+//      if further action (eg, merging, renumbering) is required.
 template<class PointList>
 label mergePoints
 (
     const PointList& points,
     const scalar mergeTol,
     const bool verbose,
-    labelList& pointMap,
-    typename PointList::const_reference origin = PointList::value_type::zero
+    labelList& pointToUnique
 );
 
 
-//- Sorts and merges points. Determines new points. Returns true if anything
-//  merged (though newPoints still sorted even if not merged).
+//- Inplace merge points, preserving the original point order.
+//- All points closer/equal mergeTol are to be merged.
+//
+//  \note The list of points must be a concrete List/Field etc.
+//  It cannot be an indirect list.
+//
+//  \param[in,out] points The unmerged points on input, and unique merged
+//      points on output
+//  \param[in] mergeTol Geometric merge tolerance
+//  \param[in] verbose Debug verbosity
+//  \param[out] pointToUnique An old-to-new mapping from original
+//      unmerged point index to the index into unique points.
+//
+//  \returns Number of point duplicates \em removed by merging.
+//  Can test as true/false to determine further actions.
 template<class PointList>
-bool mergePoints
+label inplaceMergePoints
+(
+    PointList& points,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointToUnique
+);
+
+
+//- Inplace merge points a subset of points,
+//- preserving the original point order.
+//- All points closer/equal mergeTol are to be merged.
+//
+//  \note The list of points must be a concrete List/Field etc.
+//  It cannot be an indirect list.
+//
+//  \param[in,out] points The unmerged points on input, and unique merged
+//      points on output
+//  \param[in] selection The subset of points to consider for merging
+//  \param[in] mergeTol Geometric merge tolerance
+//  \param[in] verbose Debug verbosity
+//  \param[out] pointToUnique An old-to-new mapping from original
+//      unmerged point index to the index into unique points.
+//
+//  \returns Number of point duplicates \em removed by merging.
+//  Can test as true/false to determine further actions.
+template<class PointList>
+label inplaceMergePoints
+(
+    PointList& points,
+    const labelUList& selection,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointToUnique
+);
+
+
+//- Merge points, preserving the original point order.
+//- All points closer/equal mergeTol are to be merged.
+//
+//  \param[in] points The input (unmerged) points
+//  \param[in] mergeTol Geometric merge tolerance
+//  \param[in] verbose Debug verbosity
+//  \param[out] pointToUnique An old-to-new mapping from original
+//      unmerged point index to the index into unique points.
+//  \param[out] newPoints The new unique (merged) points.
+//
+//  \returns Number of point duplicates \em removed by merging.
+//  Can test as true/false to determine further actions.
+template<class PointList>
+label mergePoints
+(
+    const PointList& points,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointToUnique,
+    List<typename PointList::value_type>& newPoints
+);
+
+
+
+//
+// Housekeeping
+//
+
+//- Deprecated (MAR-2022) reference point is ignored
+//  \deprecated(2022-03) reference point is ignored
+template<class PointList>
+FOAM_DEPRECATED_FOR(2022-03, "mergePoint() methods without reference point")
+inline label mergePoints
+(
+    const PointList& points,
+    const scalar mergeTol,
+    const bool verbose,
+    labelList& pointMap,
+    typename PointList::value_type origin  /* Used for v2112 and earlier */
+)
+{
+    return Foam::mergePoints(points, mergeTol, verbose, pointMap);
+}
+
+
+//- Deprecated (MAR-2022) reference point is ignored
+//  \deprecated(2022-03) reference point is ignored
+template<class PointList>
+FOAM_DEPRECATED_FOR(2022-03, "mergePoint() methods without reference point")
+inline label mergePoints
 (
     const PointList& points,
     const scalar mergeTol,
     const bool verbose,
     labelList& pointMap,
     List<typename PointList::value_type>& newPoints,
-    typename PointList::const_reference origin = PointList::value_type::zero
-);
+    typename PointList::value_type origin  /* Used for v2112 and earlier */
+)
+{
+    return Foam::mergePoints(points, mergeTol, verbose, pointMap, newPoints);
+}
 
 } // End namespace Foam
 
diff --git a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
index a16d392b28..c4ea5b3e3d 100644
--- a/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
+++ b/src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C
@@ -1932,20 +1932,16 @@ Foam::pointField Foam::globalMeshData::geometricSharedPoints() const
     // Merge tolerance
     scalar tolDim = matchTol_ * mesh_.bounds().mag();
 
-    // And see how many are unique
-    labelList pMap;
-    pointField mergedPoints;
-
-    Foam::mergePoints
+    labelList pointMap;
+    Foam::inplaceMergePoints
     (
         sharedPoints,   // coordinates to merge
         tolDim,         // tolerance
         false,          // verbosity
-        pMap,
-        mergedPoints
+        pointMap
     );
 
-    return mergedPoints;
+    return sharedPoints;
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
index 5272bf8069..6f6a95daf0 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2020,2021 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,11 +45,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef PatchTools_H
-#define PatchTools_H
+#ifndef Foam_PatchTools_H
+#define Foam_PatchTools_H
 
 #include "PrimitivePatch.H"
-#include "HashSet.H"
 #include "globalIndex.H"
 #include "autoPtr.H"
 
@@ -225,8 +224,12 @@ public:
     //  \param[out] mergedPoints
     //  \param[out] mergedFaces
     //  \param[out] pointMergeMap
-    //  \param[in] useLocal gather/merge the localFaces/localPoints
-    //      instead of (global) faces/points
+    //  \param[in] useLocal gather/merge patch localFaces/localPoints
+    //      instead of faces/points
+    //
+    //  \note
+    //  - OpenFOAM-v2112 and earlier: geometric merge on all patch points.
+    //  - OpenFOAM-v2206 and later: geometric merge on patch boundary points.
     template<class FaceList, class PointField>
     static void gatherAndMerge
     (
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
index 305673e989..81e587ed4b 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C
@@ -85,47 +85,71 @@ void Foam::PatchTools::gatherAndMerge
         }
     }
 
+
     // Merging points
-    bool hasMerged = false;
+    label nPointsChanged(0);
+
+    labelList boundaryPoints;
+
+    if (Pstream::parRun())
+    {
+        const globalIndex localPointAddr
+        (
+            pp.localPoints().size(),
+            globalIndex::gatherOnly{}
+        );
+
+        const globalIndex bndPointAddr
+        (
+            pp.boundaryPoints().size(),
+            globalIndex::gatherOnly{}
+        );
+
+        bndPointAddr.gather(pp.boundaryPoints(), boundaryPoints);
+
+        // Relabel according to global point offsets
+        for (const label proci : localPointAddr.subProcs())
+        {
+            SubList<label> slot(boundaryPoints, bndPointAddr.range(proci));
+            localPointAddr.inplaceToGlobal(proci, slot);
+        }
+    }
+
 
     if (Pstream::parRun() && Pstream::master())
     {
-        Field<PointType> newPoints;
-        labelList oldToNew;
+        labelList pointToUnique;
 
-        hasMerged = Foam::mergePoints
+        nPointsChanged = Foam::inplaceMergePoints
         (
             mergedPoints,
+            boundaryPoints,  // selection of points to merge
             mergeDist,
-            false,  // verbosity
-            oldToNew,
-            newPoints
+            false,           // verbose = false
+            pointToUnique
         );
 
-        if (hasMerged)
+        if (nPointsChanged)
         {
-            // Relabel faces
+            // Renumber faces to use unique point numbers
             for (auto& f : mergedFaces)
             {
-                inplaceRenumber(oldToNew, f);
+                inplaceRenumber(pointToUnique, f);
             }
 
-            // Store newly merged points
-            mergedPoints.transfer(newPoints);
-
             // Store point mapping
             if (notNull(pointMergeMap))
             {
-                pointMergeMap.transfer(oldToNew);
+                pointMergeMap.transfer(pointToUnique);
             }
         }
     }
 
-    // TDB:
-    // if (!hasMerged && notNull(pointMergeMap))
-    // {
-    //     pointMergeMap.clear();
-    // }
+    if (!nPointsChanged && notNull(pointMergeMap))
+    {
+        // Safety
+        pointMergeMap = identity(mergedPoints.size());
+    }
 }
 
 
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
index a9f47301d1..4a5d635310 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatch.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,8 +51,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef PrimitivePatch_H
-#define PrimitivePatch_H
+#ifndef Foam_PrimitivePatch_H
+#define Foam_PrimitivePatch_H
 
 #include "boolList.H"
 #include "labelList.H"
@@ -347,6 +347,8 @@ public:
             }
 
             //- Return list of boundary points, address into LOCAL point list
+            //  Uses edge addressing (if it exists) or calculates directly
+            //  from localFaces()
             const labelList& boundaryPoints() const;
 
             //- Return face-face addressing
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
index 1e9292768c..0b5fe09abc 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchAddressing.C
@@ -62,11 +62,11 @@ Foam::PrimitivePatch<FaceList, PointField>::calcAddressing() const
     // get reference to pointFaces
     const labelListList& pf = pointFaces();
 
-    // Guess the max number of edges and neighbours for a face
+    // Guess the max number of edges/neighbours for a face
     label maxEdges = 0;
     for (const auto& f : locFcs)
     {
-        maxEdges += f.size();
+        maxEdges += f.nEdges();
     }
 
     // create the lists for the various results. (resized on completion)
diff --git a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
index 5b58bcf49e..f2fdbea1c3 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/PrimitivePatch/PrimitivePatchBdryPoints.C
@@ -35,8 +35,6 @@ template<class FaceList, class PointField>
 void
 Foam::PrimitivePatch<FaceList, PointField>::calcBdryPoints() const
 {
-    DebugInFunction << "Calculating boundary points" << nl;
-
     if (boundaryPointsPtr_)
     {
         // Error to recalculate if already allocated
@@ -45,16 +43,81 @@ Foam::PrimitivePatch<FaceList, PointField>::calcBdryPoints() const
             << abort(FatalError);
     }
 
-    labelHashSet bp(2*nEdges());
+    labelHashSet bp(0);
 
-    for (const edge& e : boundaryEdges())
+    if (hasEdges())
     {
-        bp.insert(e.first());
-        bp.insert(e.second());
+        DebugInFunction
+            << "Calculating boundary points from existing addressing"
+            << nl;
+
+        bp.resize(4*nBoundaryEdges());
+
+        for (const edge& e : boundaryEdges())
+        {
+            bp.insert(e.first());
+            bp.insert(e.second());
+        }
+    }
+    else
+    {
+        DebugInFunction
+            << "Calculating boundary points with manual edge addressing"
+            << nl;
+
+
+        // Calculate manually.
+        // Needs localFaces, but uses local hashes of the edges here
+        // instead of forcing a full faceFaces/edgeFaces/faceEdges calculation
+
+        // Get reference to localFaces
+        const List<face_type>& locFcs = localFaces();
+
+        // Guess the max number of edges/neighbours for a face
+        label edgeCount = 0;
+        for (const auto& f : locFcs)
+        {
+            edgeCount += f.nEdges();
+        }
+
+        // ie, EdgeMap<label> to keep counts
+        HashTable<label, edge, Hash<edge>> knownEdges(2*edgeCount);
+
+        for (const auto& f : locFcs)
+        {
+            const label numEdges = f.nEdges();
+
+            for (label edgei = 0; edgei < numEdges; ++edgei)
+            {
+                ++ knownEdges(f.edge(edgei));
+            }
+        }
+
+        edgeCount = 0;
+
+        forAllConstIters(knownEdges, iter)
+        {
+            if (1 == iter.val())  // Singly connected edge
+            {
+                ++edgeCount;
+            }
+        }
+
+        bp.resize(4*edgeCount);
+
+        forAllConstIters(knownEdges, iter)
+        {
+            const edge& e = iter.key();
+
+            if (1 == iter.val())  // Singly connected edge
+            {
+                bp.insert(e.first());
+                bp.insert(e.second());
+            }
+        }
     }
 
     boundaryPointsPtr_.reset(new labelList(bp.sortedToc()));
-
     DebugInfo << "    Finished." << nl;
 }
 
diff --git a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
index 22acf94ab6..c03ab28114 100644
--- a/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
+++ b/src/dynamicMesh/polyMeshAdder/polyMeshAdder.C
@@ -2124,11 +2124,7 @@ Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
     //labelList sharedToMerged;
     //label nUnique = Foam::mergePoints
     //(
-    //    pointField
-    //    (
-    //        mesh.points(),
-    //        sharedPointLabels
-    //    ),
+    //    UIndirectList<point>(mesh.points(), sharedPointLabels),
     //    mergeDist,
     //    false,
     //    sharedToMerged
diff --git a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
index 74480573fe..f0c859408e 100644
--- a/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
+++ b/src/functionObjects/field/fieldValues/surfaceFieldValue/surfaceFieldValue.C
@@ -31,7 +31,6 @@ License
 #include "emptyPolyPatch.H"
 #include "coupledPolyPatch.H"
 #include "sampledSurface.H"
-#include "mergePoints.H"
 #include "indirectPrimitivePatch.H"
 #include "PatchTools.H"
 #include "addToRunTimeSelectionTable.H"
@@ -374,16 +373,25 @@ void Foam::functionObjects::fieldValues::surfaceFieldValue::combineMeshGeometry
 
     if (Pstream::parRun())
     {
-        labelList pointsMap;
-
-        PatchTools::gatherAndMerge
+        // Topological merge
+        labelList pointToGlobal;
+        labelList uniqueMeshPointLabels;
+        autoPtr<globalIndex> globalPoints;
+        autoPtr<globalIndex> globalFaces;
+        Foam::PatchTools::gatherAndMerge
         (
-            SMALL,  // mergeDist
-            pp,
-            points,
+            mesh_,
+            pp.localFaces(),
+            pp.meshPoints(),
+            pp.meshPointMap(),
+
+            pointToGlobal,
+            uniqueMeshPointLabels,
+            globalPoints,
+            globalFaces,
+
             faces,
-            pointsMap,
-            true   // useLocal=true
+            points
         );
     }
     else
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
index eef4762760..628a00400b 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C
@@ -4502,9 +4502,9 @@ void Foam::snappyLayerDriver::mapFaceZonePoints
         // too many points; the problem would be if merging baffles.
         // Trust mergeZoneBaffles to do sufficient checks.
         labelList oldToNew;
-        label nNew = mergePoints
+        label nNew = Foam::mergePoints
         (
-            pointField(mesh.points(), candidates),
+            UIndirectList<point>(mesh.points(), candidates),
             meshRefiner_.mergeDistance(),
             false,
             oldToNew
diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
index d80cc3432f..c66c3e4447 100644
--- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
+++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C
@@ -73,7 +73,7 @@ Foam::label Foam::snappySnapDriver::getCollocatedPoints
 )
 {
     labelList pointMap;
-    label nUnique = mergePoints
+    label nUnique = Foam::mergePoints
     (
         points,                         // points
         tol,                            // mergeTol
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C b/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
index 17621e0f48..565be1efea 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C
@@ -233,25 +233,22 @@ void Foam::advancingFrontAMI::distributeAndMergePatches
 
     // Merge
     labelList oldToNew;
-    pointField newTgtPoints;
-    bool hasMerged = mergePoints
+    bool nChanged = Foam::inplaceMergePoints
     (
         tgtPoints,
         SMALL,
         false,
-        oldToNew,
-        newTgtPoints
+        oldToNew
     );
 
-    if (hasMerged)
+    if (nChanged)
     {
         if (debug)
         {
-            Pout<< "Merged from " << tgtPoints.size()
-                << " down to " << newTgtPoints.size() << " points" << endl;
+            Pout<< "Merged from " << oldToNew.size()
+                << " down to " << tgtPoints.size() << " points" << endl;
         }
 
-        tgtPoints.transfer(newTgtPoints);
         for (face& f : tgtFaces)
         {
             inplaceRenumber(oldToNew, f);
diff --git a/src/meshTools/edgeMesh/edgeMesh.C b/src/meshTools/edgeMesh/edgeMesh.C
index 70b657359f..778d3eb4ce 100644
--- a/src/meshTools/edgeMesh/edgeMesh.C
+++ b/src/meshTools/edgeMesh/edgeMesh.C
@@ -212,29 +212,22 @@ void Foam::edgeMesh::scalePoints(const scalar scaleFactor)
 
 void Foam::edgeMesh::mergePoints(const scalar mergeDist)
 {
-    pointField newPoints;
     labelList pointMap;
 
-    const bool hasMerged = Foam::mergePoints
+    label nChanged = Foam::inplaceMergePoints
     (
         points_,
         mergeDist,
         false,
-        pointMap,
-        newPoints,
-        vector::zero
+        pointMap
     );
 
-    if (hasMerged)
+    if (nChanged)
     {
         pointEdgesPtr_.reset(nullptr);   // connectivity change
 
-        points_.transfer(newPoints);
-
-        forAll(edges_, edgeI)
+        for (edge& e : edges_)
         {
-            edge& e = edges_[edgeI];
-
             e[0] = pointMap[e[0]];
             e[1] = pointMap[e[1]];
         }
diff --git a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
index 6091da9d64..c2eaf27b11 100644
--- a/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
+++ b/src/meshTools/edgeMesh/extendedEdgeMesh/extendedEdgeMesh.C
@@ -1810,7 +1810,7 @@ bool Foam::extendedEdgeMesh::mergePointsAndSort
 
     // Detect and merge collocated feature points
     labelList oldToMerged;
-    label nNewPoints = ::Foam::mergePoints
+    label nNewPoints = Foam::mergePoints
     (
         points(),
         SMALL,
diff --git a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
index 9816d150c2..cb00578c98 100644
--- a/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
+++ b/src/meshTools/triSurface/booleanOps/surfaceIntersection/surfaceIntersection.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -1486,22 +1486,18 @@ const Foam::labelListList& Foam::surfaceIntersection::surf2EdgeCuts() const
 
 void Foam::surfaceIntersection::mergePoints(const scalar mergeDist)
 {
-    pointField newPoints;
     labelList pointMap;
 
-    const label nMerged = Foam::mergePoints
+    label nChanged = Foam::inplaceMergePoints
     (
         cutPoints_,
         mergeDist,
         false,
-        pointMap,
-        newPoints
+        pointMap
     );
 
-    if (nMerged)
+    if (nChanged)
     {
-        cutPoints_.transfer(newPoints);
-
         forAll(cutEdges_, edgei)
         {
             edge& e = cutEdges_[edgei];
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 2547eec57f..920ef17b67 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -1887,20 +1887,19 @@ Foam::triSurface Foam::triSurfaceTools::mergePoints
     const scalar mergeTol
 )
 {
-    pointField newPoints(surf.nPoints());
+    labelList pointMap;
+    labelList uniquePoints;
 
-    labelList pointMap(surf.nPoints());
-
-    bool hasMerged = Foam::mergePoints
+    label nChanged = Foam::mergePoints
     (
         surf.localPoints(),
-        mergeTol,
-        false,
         pointMap,
-        newPoints
+        uniquePoints,
+        mergeTol,
+        false
     );
 
-    if (hasMerged)
+    if (nChanged)
     {
         // Pack the triangles
 
@@ -1923,6 +1922,8 @@ Foam::triSurface Foam::triSurfaceTools::mergePoints
         }
         newTriangles.resize(nNewTris);
 
+        pointField newPoints(surf.localPoints(), uniquePoints);
+
         return triSurface
         (
             newTriangles,
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 6fd356690f..b7285cc19a 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -1860,9 +1860,8 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 //            (
 //                globalBorderCentres,
 //                mergeDist_,
-//                false,          //const bool verbose,
+//                false,          // verbose = false
 //                allToMerged
-//                // maybe bounds().mid() ?
 //            );
 //
 //            if (debug)
@@ -2006,9 +2005,8 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
                 (
                     allCentres,
                     mergeDist_,
-                    false,          //const bool verbose,
+                    false,          // verbose = false
                     allToMerged
-                    // maybe bounds().mid() ?
                 );
 
                 if (debug)
diff --git a/src/sampling/meshToMesh/meshToMeshParallelOps.C b/src/sampling/meshToMesh/meshToMeshParallelOps.C
index 6f47496b5d..4f15b95ba7 100644
--- a/src/sampling/meshToMesh/meshToMeshParallelOps.C
+++ b/src/sampling/meshToMesh/meshToMeshParallelOps.C
@@ -862,25 +862,22 @@ void Foam::meshToMesh::distributeAndMergeCells
         // only merging points in debug mode
 
         labelList oldToNew;
-        pointField newTgtPoints;
-        bool hasMerged = mergePoints
+        label nChanged = Foam::inplaceMergePoints
         (
             tgtPoints,
             SMALL,
             false,
-            oldToNew,
-            newTgtPoints
+            oldToNew
         );
 
-        if (hasMerged)
+        if (nChanged)
         {
-            Pout<< "Merged from " << tgtPoints.size()
-                << " down to " << newTgtPoints.size() << " points" << endl;
+            Pout<< "Merged from " << oldToNew.size()
+                << " down to " << tgtPoints.size() << " points" << endl;
 
-            tgtPoints.transfer(newTgtPoints);
-            forAll(tgtFaces, i)
+            for (auto& f : tgtFaces)
             {
-                inplaceRenumber(oldToNew, tgtFaces[i]);
+                inplaceRenumber(oldToNew, f);
             }
         }
     }
diff --git a/src/sampling/sampledSet/patchEdge/patchEdgeSet.C b/src/sampling/sampledSet/patchEdge/patchEdgeSet.C
index 3475a03049..173f78b5bf 100644
--- a/src/sampling/sampledSet/patchEdge/patchEdgeSet.C
+++ b/src/sampling/sampledSet/patchEdge/patchEdgeSet.C
@@ -131,13 +131,12 @@ void Foam::patchEdgeSet::genSamples()
 
 
     labelList pointMap;
-    const label nMerged = mergePoints
+    const label nMerged = Foam::mergePoints
     (
         samplingPts,
-        SMALL,          //const scalar mergeTol,
-        false,          //const bool verbose,
-        pointMap,
-        origin_
+        SMALL,          // mergeTol
+        false,          // verbose = false
+        pointMap
     );
 
     if (nMerged == samplingPts.size())
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.C b/src/sampling/surface/isoSurface/isoSurfaceCell.C
index 434e348c43..5c94e349ae 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceCell.C
@@ -767,7 +767,7 @@ Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
     }
 
     pointField newPoints;
-    mergePoints
+    Foam::mergePoints
     (
         triPoints,
         mergeDistance_,
@@ -779,27 +779,25 @@ Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
     // Check that enough merged.
     if (debug)
     {
-        Pout<< "isoSurfaceCell : merged from " << triPoints.size()
+        Pout<< "isoSurfaceCell : merged from " << triPointReverseMap.size()
             << " points down to " << newPoints.size() << endl;
 
-        pointField newNewPoints;
         labelList oldToNew;
-        bool hasMerged = mergePoints
+        const label nUnique = Foam::mergePoints
         (
             newPoints,
             mergeDistance_,
             true,
-            oldToNew,
-            newNewPoints
+            oldToNew
         );
 
-        if (hasMerged)
+        if (nUnique != newPoints.size())
         {
             FatalErrorInFunction
                 << "Merged points contain duplicates"
                 << " when merging with distance " << mergeDistance_ << endl
                 << "merged:" << newPoints.size() << " re-merged:"
-                << newNewPoints.size()
+                << nUnique
                 << abort(FatalError);
         }
     }
@@ -852,30 +850,28 @@ Foam::triSurface Foam::isoSurfaceCell::stitchTriPoints
             centres[triI] = tris[triI].centre(newPoints);
         }
 
-        pointField mergedCentres;
         labelList oldToMerged;
-        bool hasMerged = mergePoints
+        label nUnique = Foam::mergePoints
         (
             centres,
             mergeDistance_,
             false,
-            oldToMerged,
-            mergedCentres
+            oldToMerged
         );
 
         if (debug)
         {
             Pout<< "isoSurfaceCell : detected "
-                << centres.size()-mergedCentres.size()
+                << (oldToMerged.size() - nUnique)
                 << " duplicate triangles." << endl;
         }
 
-        if (hasMerged)
+        if (oldToMerged.size() != nUnique)
         {
             // Filter out duplicates.
             label newTriI = 0;
             DynamicList<label> newToOldTri(tris.size());
-            labelList newToMaster(mergedCentres.size(), -1);
+            labelList newToMaster(nUnique, -1);
             forAll(tris, triI)
             {
                 label mergedI = oldToMerged[triI];
@@ -918,26 +914,25 @@ void Foam::isoSurfaceCell::calcAddressing
         edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
     }
 
-    pointField mergedCentres;
     labelList oldToMerged;
-    bool hasMerged = mergePoints
+    const label nUnique = Foam::mergePoints
     (
         edgeCentres,
         mergeDistance_,
         false,
-        oldToMerged,
-        mergedCentres
+        oldToMerged
     );
 
     if (debug)
     {
         Pout<< "isoSurfaceCell : detected "
-            << mergedCentres.size()
+            << nUnique
             << " edges on " << surf.size() << " triangles." << endl;
     }
 
-    if (!hasMerged)
+    if (nUnique == edgeCentres.size())
     {
+        // Nothing to do
         return;
     }
 
@@ -954,9 +949,9 @@ void Foam::isoSurfaceCell::calcAddressing
 
 
     // Determine edgeFaces
-    edgeFace0.setSize(mergedCentres.size());
+    edgeFace0.resize(nUnique);
     edgeFace0 = -1;
-    edgeFace1.setSize(mergedCentres.size());
+    edgeFace1.resize(nUnique);
     edgeFace1 = -1;
     edgeFacesRest.clear();
 
@@ -976,7 +971,7 @@ void Foam::isoSurfaceCell::calcAddressing
         else
         {
             //WarningInFunction
-            //    << "Edge " << edgeI << " with centre " << mergedCentres[edgeI]
+            //    << "Edge " << edgeI << " with centre "
             //    << " used by more than two triangles: " << edgeFace0[edgeI]
             //    << ", "
             //    << edgeFace1[edgeI] << " and " << triI << endl;
diff --git a/src/sampling/surface/isoSurface/isoSurfacePoint.C b/src/sampling/surface/isoSurface/isoSurfacePoint.C
index b0d3f6be40..7d345270bb 100644
--- a/src/sampling/surface/isoSurface/isoSurfacePoint.C
+++ b/src/sampling/surface/isoSurface/isoSurfacePoint.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2021 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -906,7 +906,7 @@ Foam::triSurface Foam::isoSurfacePoint::stitchTriPoints
     // Check that enough merged.
     if (debug)
     {
-        Pout<< "isoSurfacePoint : merged from " << triPoints.size()
+        Pout<< "isoSurfacePoint : merged from " << triPointReverseMap.size()
             << " down to " << newPoints.size() << " unique points." << endl;
 
         pointField newNewPoints;
@@ -925,7 +925,7 @@ Foam::triSurface Foam::isoSurfacePoint::stitchTriPoints
             FatalErrorInFunction
                 << "Merged points contain duplicates"
                 << " when merging with distance " << mergeDistance_ << endl
-                << "merged:" << newPoints.size() << " re-merged:"
+                << "merged:" << oldToNew.size() << " re-merged:"
                 << newNewPoints.size()
                 << abort(FatalError);
         }
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index ef90d67771..c770aed882 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -695,15 +695,13 @@ bool Foam::MeshedSurface<Face>::stitchFaces
     const bool verbose
 )
 {
-    pointField& pointLst = this->storedPoints();
+    pointField& ps = this->storedPoints();
 
-    // Merge points
-    labelList  pointMap(pointLst.size());
-    pointField newPoints(pointLst.size());
+    // Merge points (inplace)
+    labelList pointMap;
+    label nChanged = Foam::inplaceMergePoints(ps, tol, verbose, pointMap);
 
-    bool hasMerged = mergePoints(pointLst, tol, verbose, pointMap, newPoints);
-
-    if (!hasMerged)
+    if (!nChanged)
     {
         return false;
     }
@@ -713,9 +711,6 @@ bool Foam::MeshedSurface<Face>::stitchFaces
         InfoInFunction<< "Renumbering all faces" << endl;
     }
 
-    // Set the coordinates to the merged ones
-    pointLst.transfer(newPoints);
-
     List<Face>& faceLst = this->storedFaces();
 
     labelList faceMap(faceLst.size(), -1);
diff --git a/src/surfMesh/surfaceFormats/tri/TRIReader.C b/src/surfMesh/surfaceFormats/tri/TRIReader.C
index 8125cc497d..d1e42d6e1c 100644
--- a/src/surfMesh/surfaceFormats/tri/TRIReader.C
+++ b/src/surfMesh/surfaceFormats/tri/TRIReader.C
@@ -31,6 +31,7 @@ License
 #include "IFstream.H"
 #include "IOmanip.H"
 #include "StringStream.H"
+#include "DynamicList.H"
 #include "mergePoints.H"
 #include "Map.H"
 
diff --git a/src/surfMesh/triSurface/triSurfaceStitch.C b/src/surfMesh/triSurface/triSurfaceStitch.C
index e1fcb468de..7d179da1f4 100644
--- a/src/surfMesh/triSurface/triSurfaceStitch.C
+++ b/src/surfMesh/triSurface/triSurfaceStitch.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,24 +38,20 @@ bool Foam::triSurface::stitchTriangles
     bool verbose
 )
 {
-    pointField& ps = storedPoints();
+    pointField& ps = this->storedPoints();
 
-    // Merge points
+    // Merge points (inplace)
     labelList pointMap;
-    pointField newPoints;
-    bool hasMerged = mergePoints(ps, tol, verbose, pointMap, newPoints);
+    label nChanged = Foam::inplaceMergePoints(ps, tol, verbose, pointMap);
 
-    if (hasMerged)
+    if (nChanged)
     {
         if (verbose)
         {
-            Pout<< "stitchTriangles : Merged from " << ps.size()
-                << " points down to " << newPoints.size() << endl;
+            Pout<< "stitchTriangles : Merged from " << pointMap.size()
+                << " points down to " << ps.size() << endl;
         }
 
-        // Set the coordinates to the merged ones
-        ps.transfer(newPoints);
-
         // Reset the triangle point labels to the unique points array
         label newTriangleI = 0;
         forAll(*this, i)
@@ -143,7 +140,7 @@ bool Foam::triSurface::stitchTriangles
         }
     }
 
-    return hasMerged;
+    return bool(nChanged);
 }