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);
 }