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

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,

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 "ListOps.H"  // sortedOrder, ListOps::identity
-#include "point.H"
-#include "Field.H"
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Implementation  * * * * * * * * * * * * * * //
 
-template<class PointList>
+template<class PointList, class IndexerOp>
-Foam::label Foam::mergePoints
+Foam::label Foam::Detail::mergePoints
 (
     const PointList& points,
+    const IndexerOp& indexer,
+    const label nSubPoints,
+    labelList& pointToUnique,
+    labelList& uniquePoints,
     const scalar mergeTol,
-    const bool verbose,
+    const bool verbose
-    labelList& pointMap,
-    typename PointList::const_reference origin
 )
 {
-    typedef typename PointList::value_type point_type;
+    const label nTotPoints = points.size();
 
-    const label nPoints = points.size();
+    if (!nTotPoints || !nSubPoints)
-
-    // Create an old to new point mapping array
-    pointMap.resize_nocopy(nPoints);
-    pointMap = -1;
-
-    if (!nPoints)
     {
-        return 0;
+        // Nothing to do
+        pointToUnique = identity(nTotPoints);
+        uniquePoints = pointToUnique;
+        return 0;  // No points removed
     }
 
-    point_type compareOrigin = origin;
+    // Properly size for old to new mapping array
-    if (origin == point_type::max)
+    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.
+        comparePoint = min(comparePoint, points[indexer(pointi)]);
-        // 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;
     }
 
-    // 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
+    // Use magSqr distance for the points sort order
-    List<scalar> magSqrDist(nPoints);
+    List<scalar> sqrDist(nSubPoints);
-    forAll(points, pointi)
+    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));
 
-
+    List<scalar> sortedTol(nSubPoints);
-    Field<scalar> sortedTol(nPoints);
     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()))
+                // Use scalar precision
-              + mag(scalar(pt.y() - compareOrigin.y()))
+                mag(scalar(p.x() - comparePoint.x()))
-              + mag(scalar(pt.z() - compareOrigin.z()))
+              + mag(scalar(p.y() - comparePoint.y()))
+              + mag(scalar(p.z() - comparePoint.z()))
+            )
         );
     }
 
-    label newPointi = 0;
 
-    // Handle 0th point separately (is always unique)
+    // Bookkeeping parameters
-    label pointi = order[0];
+    // ~~~~~~~~~~~~~~~~~~~~~~
-    pointMap[pointi] = newPointi++;
 
-    /// if (verbose)
+    // Will only be working on a subset of the points
-    /// {
+    // Can use a slice of pointToUnique (full length not needed until later).
-    ///     Pout<< "Foam::mergePoints : [0] Uniq point " << pointi << endl;
+    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
+        // Compare to previous points to find equal one
-        // NOTE: not yet using point_type template parameter
+        // - automatically a no-op for sorti == 0 (the first point)
-        const point pt
-        (
-            scalar(points[pointi].x()),
-            scalar(points[pointi].y()),
-            scalar(points[pointi].z())
-        );
 
-
+        bool matched = false;
-        // Compare to previous points to find equal one.
-        label equalPointi = -1;
 
         for
         (
             label prevSorti = sorti - 1;
+            (
                 prevSorti >= 0
-         && (mag(magSqrDist[order[prevSorti]] - mag2) <= sortedTol[sorti]);
+             && (mag(sqrDist[order[prevSorti]] - currDist) <= sortedTol[sorti])
+            );
             --prevSorti
         )
         {
             const label prevPointi = order[prevSorti];
+            const auto& prevPoint = points[indexer(prevPointi)];
 
-            // Convert to scalar precision
+            // Matched within tolerance?
-            // NOTE: not yet using point_type template parameter
+            matched =
-            const point prevPt
             (
-                scalar(points[prevPointi].x()),
+                (
-                scalar(points[prevPointi].y()),
+                    // Use scalar precision
-                scalar(points[prevPointi].z())
+                    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.
+                // Both pointi and prevPointi have similar coordinates.
-                equalPointi = prevPointi;
+                // 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 (!matched)
-        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
         {
             // Differs. Store new point.
+            subPointMap[pointi] = nNewPoints++;
 
+            /// Too verbose
             ///if (verbose)
             ///{
-            ///     Pout<< "Foam::mergePoints : [" << newPointi
+            ///    Pout<< "Foam::mergePoints : [" << subPointMap[pointi]
-            ///         << "] Uniq point " << pointi << endl;
+            ///        << "] Point " << pointi << endl;
             ///}
+        }
+        ++newPointCounts[subPointMap[pointi]];
+    }
 
-            pointMap[pointi] = newPointi++;
+    const label nDupPoints(nSubPoints - nNewPoints);
-        }
+    const label nUniqPoints(nTotPoints - nDupPoints);
-    }
 
     if (verbose)
     {
         Pout<< "Foam::mergePoints : "
-            << newPointi << " of " << points.size() << " unique points"
+            << "Merging removed " << nDupPoints << '/'
-            << endl;
+            << 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,
+    labelList& pointToUnique
-    List<typename PointList::value_type>& newPoints,
-    typename PointList::const_reference origin
 )
 {
-    const label nUnique = Foam::mergePoints
+    labelList uniquePoints;
+    const label nChanged = Foam::mergePoints
     (
         points,
+        pointToUnique,
+        uniquePoints,
         mergeTol,
-        verbose,
+        verbose
-        pointMap,
-        origin
     );
 
-    newPoints.resize_nocopy(nUnique);
+    // Number of unique points
-    forAll(pointMap, pointi)
+    return (points.size() - nChanged);
-    {
-        newPoints[pointMap[pointi]] = points[pointi];
 }
 
-    return (nUnique != points.size());
+
+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)
+    {
+        // Overwrite
+        points = List<typename PointList::value_type>(points, uniquePoints);
+    }
+    else
+    {
+        // TDB:
+        // pointToUnique.clear();
+    }
+
+    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;
 }
 
 

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
+#ifndef Foam_mergePoints_H
-#define 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.
+//- Implementation detail for Foam::mergePoints
-//  Returns the number of unique points and a map from old to new.
+//
+//  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,
+    labelList& pointToUnique
-    typename PointList::const_reference origin = PointList::value_type::zero
 );
 
 
-//- Sorts and merges points. Determines new points. Returns true if anything
+//- Inplace merge points, preserving the original point order.
-//  merged (though newPoints still sorted even if not merged).
+//- 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
 

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 pointMap;
-    labelList pMap;
+    Foam::inplaceMergePoints
-    pointField mergedPoints;
-
-    Foam::mergePoints
     (
         sharedPoints,   // coordinates to merge
         tolDim,         // tolerance
         false,          // verbosity
-        pMap,
+        pointMap
-        mergedPoints
     );
 
-    return mergedPoints;
+    return sharedPoints;
 }
 
 

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
+#ifndef Foam_PatchTools_H
-#define 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
+    //  \param[in] useLocal gather/merge patch localFaces/localPoints
-    //      instead of (global) faces/points
+    //      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
     (

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 pointToUnique;
-        labelList oldToNew;
 
-        hasMerged = Foam::mergePoints
+        nPointsChanged = Foam::inplaceMergePoints
         (
             mergedPoints,
+            boundaryPoints,  // selection of points to merge
             mergeDist,
-            false,  // verbosity
+            false,           // verbose = false
-            oldToNew,
+            pointToUnique
-            newPoints
         );
 
-        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 (!nPointsChanged && notNull(pointMergeMap))
-    // if (!hasMerged && notNull(pointMergeMap))
+    {
-    // {
+        // Safety
-    //     pointMergeMap.clear();
+        pointMergeMap = identity(mergedPoints.size());
-    // }
+    }
 }
 
 

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
+#ifndef Foam_PrimitivePatch_H
-#define 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

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)

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);
+
+    if (hasEdges())
+    {
+        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;
 }
 

src/dynamicMesh/polyMeshAdder/polyMeshAdder.C

@@ -2124,11 +2124,7 @@ Foam::Map<Foam::label> Foam::polyMeshAdder::findSharedPoints
     //labelList sharedToMerged;
     //label nUnique = Foam::mergePoints
     //(
-    //    pointField
+    //    UIndirectList<point>(mesh.points(), sharedPointLabels),
-    //    (
-    //        mesh.points(),
-    //        sharedPointLabels
-    //    ),
     //    mergeDist,
     //    false,
     //    sharedToMerged

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;
+        // Topological merge
-
+        labelList pointToGlobal;
-        PatchTools::gatherAndMerge
+        labelList uniqueMeshPointLabels;
+        autoPtr<globalIndex> globalPoints;
+        autoPtr<globalIndex> globalFaces;
+        Foam::PatchTools::gatherAndMerge
         (
-            SMALL,  // mergeDist
+            mesh_,
-            pp,
+            pp.localFaces(),
-            points,
+            pp.meshPoints(),
+            pp.meshPointMap(),
+
+            pointToGlobal,
+            uniqueMeshPointLabels,
+            globalPoints,
+            globalFaces,
+
             faces,
-            pointsMap,
+            points
-            true   // useLocal=true
         );
     }
     else

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

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

src/meshTools/AMIInterpolation/AMIInterpolation/advancingFrontAMI/advancingFrontAMIParallelOps.C

@@ -233,25 +233,22 @@ void Foam::advancingFrontAMI::distributeAndMergePatches
 
     // Merge
     labelList oldToNew;
-    pointField newTgtPoints;
+    bool nChanged = Foam::inplaceMergePoints
-    bool hasMerged = mergePoints
     (
         tgtPoints,
         SMALL,
         false,
-        oldToNew,
+        oldToNew
-        newTgtPoints
     );
 
-    if (hasMerged)
+    if (nChanged)
     {
         if (debug)
         {
-            Pout<< "Merged from " << tgtPoints.size()
+            Pout<< "Merged from " << oldToNew.size()
-                << " down to " << newTgtPoints.size() << " points" << endl;
+                << " down to " << tgtPoints.size() << " points" << endl;
         }
 
-        tgtPoints.transfer(newTgtPoints);
         for (face& f : tgtFaces)
         {
             inplaceRenumber(oldToNew, f);

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,
+        pointMap
-        newPoints,
-        vector::zero
     );
 
-    if (hasMerged)
+    if (nChanged)
     {
         pointEdgesPtr_.reset(nullptr);   // connectivity change
 
-        points_.transfer(newPoints);
+        for (edge& e : edges_)
-
-        forAll(edges_, edgeI)
         {
-            edge& e = edges_[edgeI];
-
             e[0] = pointMap[e[0]];
             e[1] = pointMap[e[1]];
         }

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,

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,
+        pointMap
-        newPoints
     );
 
-    if (nMerged)
+    if (nChanged)
     {
-        cutPoints_.transfer(newPoints);
-
         forAll(cutEdges_, edgei)
         {
             edge& e = cutEdges_[edgei];

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());
+    label nChanged = Foam::mergePoints
-
-    bool hasMerged = 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,

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)

src/sampling/meshToMesh/meshToMeshParallelOps.C

@@ -862,25 +862,22 @@ void Foam::meshToMesh::distributeAndMergeCells
         // only merging points in debug mode
 
         labelList oldToNew;
-        pointField newTgtPoints;
+        label nChanged = Foam::inplaceMergePoints
-        bool hasMerged = mergePoints
         (
             tgtPoints,
             SMALL,
             false,
-            oldToNew,
+            oldToNew
-            newTgtPoints
         );
 
-        if (hasMerged)
+        if (nChanged)
         {
-            Pout<< "Merged from " << tgtPoints.size()
+            Pout<< "Merged from " << oldToNew.size()
-                << " down to " << newTgtPoints.size() << " points" << endl;
+                << " down to " << tgtPoints.size() << " points" << endl;
 
-            tgtPoints.transfer(newTgtPoints);
+            for (auto& f : tgtFaces)
-            forAll(tgtFaces, i)
             {
-                inplaceRenumber(oldToNew, tgtFaces[i]);
+                inplaceRenumber(oldToNew, f);
             }
         }
     }

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,
+        SMALL,          // mergeTol
-        false,          //const bool verbose,
+        false,          // verbose = false
-        pointMap,
+        pointMap
-        origin_
     );
 
     if (nMerged == samplingPts.size())

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,
+            oldToNew
-            newNewPoints
         );
 
-        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,
+            oldToMerged
-            mergedCentres
         );
 
         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,
+        oldToMerged
-        mergedCentres
     );
 
     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;

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

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
+    // Merge points (inplace)
-    labelList  pointMap(pointLst.size());
+    labelList pointMap;
-    pointField newPoints(pointLst.size());
+    label nChanged = Foam::inplaceMergePoints(ps, tol, verbose, pointMap);
 
-    bool hasMerged = mergePoints(pointLst, tol, verbose, pointMap, newPoints);
+    if (!nChanged)
-
-    if (!hasMerged)
     {
         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);

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"
 

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;
+    label nChanged = Foam::inplaceMergePoints(ps, tol, verbose, pointMap);
-    bool hasMerged = mergePoints(ps, tol, verbose, pointMap, newPoints);
 
-    if (hasMerged)
+    if (nChanged)
     {
         if (verbose)
         {
-            Pout<< "stitchTriangles : Merged from " << ps.size()
+            Pout<< "stitchTriangles : Merged from " << pointMap.size()
-                << " points down to " << newPoints.size() << endl;
+                << " 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);
 }