diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index 0d7f9ec25a17876426aedf50bb0f3a81ab60484c..37433f884a013ac87620804b18d8a376c6895a52 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-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -278,7 +278,9 @@ Foam::MeshedSurface<Face>::MeshedSurface
     MeshReference(faceLst, pointLst), // Copy construct
     faceIds_(),
     zones_(zoneLst)
-{}
+{
+    this->checkZones(false);  // Non-verbose fix zones
+}
 
 
 template<class Face>
@@ -292,7 +294,9 @@ Foam::MeshedSurface<Face>::MeshedSurface
     MeshReference(faceLst, pointLst, true), // Move construct
     faceIds_(),
     zones_(zoneLst)
-{}
+{
+    this->checkZones(false);  // Non-verbose fix zones
+}
 
 
 template<class Face>
@@ -607,7 +611,7 @@ void Foam::MeshedSurface<Face>::movePoints(const pointField& newPoints)
 {
     MeshReference::clearGeom();  // Changes areas, normals etc.
 
-    // Adapt for new point position
+    // Adapt for new point positions
     MeshReference::movePoints(newPoints);
 
     // Copy new points
@@ -618,17 +622,16 @@ void Foam::MeshedSurface<Face>::movePoints(const pointField& newPoints)
 template<class Face>
 void Foam::MeshedSurface<Face>::scalePoints(const scalar scaleFactor)
 {
-    // Avoid bad scaling
-    if (scaleFactor > 0 && scaleFactor != 1.0)
+    // Avoid bad or no scaling
+    if (scaleFactor > SMALL && !equal(scaleFactor, 1))
     {
-        MeshReference::clearGeom();  // Changes areas, normals etc.
+        // Remove all geometry dependent data
+        this->clearTopology();
 
-        pointField newPoints(scaleFactor*this->points());
+        // Adapt for new point positions
+        MeshReference::movePoints(pointField());
 
-        // Adapt for new point position
-        MeshReference::movePoints(newPoints);
-
-        storedPoints() = std::move(newPoints);
+        this->storedPoints() *= scaleFactor;
     }
 }
 
@@ -645,6 +648,46 @@ void Foam::MeshedSurface<Face>::cleanup(const bool verbose)
 }
 
 
+template<class Face>
+void Foam::MeshedSurface<Face>::compactPoints(labelList& pointMap)
+{
+    this->clearOut();   // Topology changes
+
+    // Remove unused points while walking and renumbering faces
+    // in visit order - walk order as per localFaces()
+
+    labelList oldToCompact(this->points().size(), -1);
+    DynamicList<label> compactPointMap(oldToCompact.size());
+
+    for (auto& f : this->storedFaces())
+    {
+        for (label& pointi : f)
+        {
+            label compacti = oldToCompact[pointi];
+            if (compacti == -1)
+            {
+                compacti = compactPointMap.size();
+                oldToCompact[pointi] = compacti;
+                compactPointMap.append(pointi);
+            }
+            pointi = compacti;
+        }
+    }
+
+    pointField newPoints
+    (
+        UIndirectList<point>(this->points(), compactPointMap)
+    );
+
+    this->swapPoints(newPoints);
+
+    if (notNull(pointMap))
+    {
+        pointMap.transfer(compactPointMap);
+    }
+}
+
+
 template<class Face>
 bool Foam::MeshedSurface<Face>::stitchFaces
 (
@@ -1345,7 +1388,7 @@ void Foam::MeshedSurface<Face>::transfer
     (
         std::move(surf.storedPoints()),
         std::move(faceLst),
-        std::move(zoneLst)
+        zoneLst
     );
 
     surf.clear();
@@ -1370,13 +1413,16 @@ void Foam::MeshedSurface<Face>::swapFaces(List<Face>& faces)
     this->storedFaceIds().clear();  // Likely to be invalid
 
     this->storedFaces().swap(faces);
+
+    this->checkZones(false);  // Non-verbose fix zones
 }
 
 
 template<class Face>
 void Foam::MeshedSurface<Face>::swapPoints(pointField& points)
 {
-    MeshReference::clearOut();  // Topology changes
+    // Adapt for new point positions
+    MeshReference::movePoints(points);
 
     this->storedPoints().swap(points);
 }
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.H b/src/surfMesh/MeshedSurface/MeshedSurface.H
index 5da16acea27e348f81f6ad57cb79479f89aab732..77f7da7c8f867cfdbbc82d8443647de9047bd8f2 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.H
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -159,8 +159,10 @@ protected:
         //  No general form, only specializations.
         void transcribe(MeshedSurface<face>& surf);
 
-        //- Basic sanity check on zones
-        void checkZones();
+        //- Sanity check/resizing on zones.
+        //  Adjust zones so that they cover the number of faces
+        //  The last zone will be extended as needed
+        void checkZones(const bool verbose = true);
 
         //- Non-const access to global points
         pointField& storedPoints()
@@ -496,6 +498,14 @@ public:
         //- Remove invalid faces
         virtual void cleanup(const bool verbose);
 
+        //- Remove unused points and renumber faces in local visit order
+        //
+        //  \param[out] pointMap from new to old points (optional)
+        virtual void compactPoints
+        (
+            labelList& pointMap = const_cast<labelList&>(labelList::null())
+        );
+
         virtual bool stitchFaces
         (
             const scalar tol=SMALL,
diff --git a/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C b/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
index 051a413cbec0ec093f8f696d31012d4bdf115f1b..116771bcb7cf3c563cddce461cbaf372c711cad5 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurfaceZones.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,35 +32,55 @@ License
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class Face>
-void Foam::MeshedSurface<Face>::checkZones()
+void Foam::MeshedSurface<Face>::checkZones(const bool verbose)
 {
-    // extra safety, ensure we have at some zones
-    // and they cover all the faces - fix start silently
-
     auto& zones = this->storedZones();
 
-    label count = 0;
+    // Check that zones (if any) they cover the faces
+    // - fix the start silently
+
+    bool zonesTooBig(false);
+
+    const label maxCount = this->size();
+
+    label start = 0;
     for (surfZone& zn : zones)
     {
-        zn.start() = count;
-        count += zn.size();
+        zn.start() = start;
+        start += zn.size();
+        if (start > maxCount)
+        {
+            zonesTooBig = true;  // Zones exceed what is required
+            zn.size() = (maxCount - zn.start());
+            start = (zn.start() + zn.size());
+        }
     }
 
     if (!zones.empty())
     {
-        if (count < this->size())
-        {
-            WarningInFunction
-                << "more faces " << this->size() << " than zones " << count
-                << " ... extending final zone" << nl;
+        surfZone& zn = zones.last();
 
-            zones.last().size() += count - this->size();
+        if ((zn.start() + zn.size()) < maxCount)
+        {
+            // Zones address less than expected - extend final zone
+            zn.size() += maxCount - zn.start();
+
+            if (verbose)
+            {
+                WarningInFunction
+                    << "Surface has more faces " << maxCount
+                    << " than zone addressing ... extending final zone" << nl;
+            }
         }
-        else if (count > this->size())
+        else if (zonesTooBig)
         {
-            FatalErrorInFunction
-                << "more zones " << count << " than faces " << this->size()
-                << exit(FatalError);
+            if (verbose)
+            {
+                WarningInFunction
+                    << "Surface has more zone addressing than faces "
+                    << maxCount
+                    << " ... trucated/resized accordingly" << nl;
+            }
         }
     }
 }
diff --git a/src/surfMesh/ModifiableMeshedSurface/ModifiableMeshedSurface.H b/src/surfMesh/ModifiableMeshedSurface/ModifiableMeshedSurface.H
index 35516883fa3d3656d6e7e1f70c9151dd186834a1..41f4cd70e60368767ec91107cc2d30853eb303aa 100644
--- a/src/surfMesh/ModifiableMeshedSurface/ModifiableMeshedSurface.H
+++ b/src/surfMesh/ModifiableMeshedSurface/ModifiableMeshedSurface.H
@@ -58,7 +58,7 @@ public:
 
     // Constructors
 
-        //- Construct null. Use swap or transfer to manage content
+        //- Default construct. Use swap or transfer to manage content
         ModifiableMeshedSurface()
         :
             MeshedSurface<Face>()
diff --git a/src/surfMesh/polySurface/polySurface.H b/src/surfMesh/polySurface/polySurface.H
index 98fa34f72563db01a99c2f3815c22f1964178172..9fb5479b8011e598245098a74ba563da5c87653d 100644
--- a/src/surfMesh/polySurface/polySurface.H
+++ b/src/surfMesh/polySurface/polySurface.H
@@ -249,19 +249,19 @@ public:
         }
 
         //- Return face area vectors (normals)
-        inline const vectorField& Sf() const
+        const vectorField& Sf() const
         {
             return MeshReference::faceAreas();
         }
 
         //- Return face area magnitudes
-        inline const scalarField& magSf() const
+        const scalarField& magSf() const
         {
             return MeshReference::magFaceAreas();
         }
 
         //- Face centres
-        inline const vectorField& Cf() const
+        const vectorField& Cf() const
         {
             return MeshReference::faceCentres();
         }
diff --git a/src/surfMesh/surfMesh/surfMesh.C b/src/surfMesh/surfMesh/surfMesh.C
index e41911fefe591495a25822f3a1be79a0db60a462..6daed2fbc9dd288e62496dd3e90584760a29e315 100644
--- a/src/surfMesh/surfMesh/surfMesh.C
+++ b/src/surfMesh/surfMesh/surfMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -465,38 +465,62 @@ const Foam::faceList& Foam::surfMesh::faces() const
 }
 
 
-void Foam::surfMesh::checkZones()
+void Foam::surfMesh::checkZones(const bool verbose)
 {
-    // Extra safety, ensure we have at some zones
-    // and they cover all the faces - fix start silently
+    auto& zones = this->storedZones();
 
-    if (surfZones_.size() <= 1)
+    if (zones.size() <= 1)
     {
         removeZones();
         return;
     }
 
-    label count = 0;
-    for (surfZone& zn : surfZones_)
-    {
-        zn.start() = count;
-        count += zn.size();
-    }
+    // Check that zones cover all faces exactly
+    // - fix each start silently
 
-    if (count < nFaces())
-    {
-        WarningInFunction
-            << "More faces " << nFaces() << " than zones " << count
-            << " ... extending final zone"
-            << endl;
+    bool zonesTooBig(false);
 
-        surfZones_.last().size() += count - nFaces();
+    const label maxCount = this->nFaces();
+
+    label start = 0;
+    for (surfZone& zn : zones)
+    {
+        zn.start() = start;
+        start += zn.size();
+        if (start > maxCount)
+        {
+            zonesTooBig = true;  // Zones exceed what is required
+            zn.size() = (maxCount - zn.start());
+            start = (zn.start() + zn.size());
+        }
     }
-    else if (size() < count)
+
+    if (!zones.empty())
     {
-        FatalErrorInFunction
-            << "More zones " << count << " than faces " << nFaces()
-            << exit(FatalError);
+        surfZone& zn = zones.last();
+
+        if ((zn.start() + zn.size()) < maxCount)
+        {
+            // Zones address less than expected - extend final zone
+            zn.size() += maxCount - zn.start();
+
+            if (verbose)
+            {
+                WarningInFunction
+                    << "Surface has more faces " << maxCount
+                    << " than zone addressing ... extending final zone" << nl;
+            }
+        }
+        else if (zonesTooBig)
+        {
+            if (verbose)
+            {
+                WarningInFunction
+                    << "Surface has more zone addressing than faces "
+                    << maxCount
+                    << " ... trucated/resized accordingly" << nl;
+            }
+        }
     }
 }
 
diff --git a/src/surfMesh/surfMesh/surfMesh.H b/src/surfMesh/surfMesh/surfMesh.H
index 70e0bad1e0398fce37b2e09d43e549ed660ecbd5..9994e445d56f25bfa57f8dbd704c230a5b5d99ab 100644
--- a/src/surfMesh/surfMesh/surfMesh.H
+++ b/src/surfMesh/surfMesh/surfMesh.H
@@ -241,19 +241,19 @@ public:
         }
 
         //- Return face area vectors (normals)
-        inline const vectorField& Sf() const
+        const vectorField& Sf() const
         {
             return MeshReference::faceAreas();
         }
 
         //- Return face area magnitudes
-        inline const scalarField& magSf() const
+        const scalarField& magSf() const
         {
             return MeshReference::magFaceAreas();
         }
 
         //- Face centres
-        inline const vectorField& Cf() const
+        const vectorField& Cf() const
         {
             return MeshReference::faceCentres();
         }
@@ -272,7 +272,7 @@ public:
         void removeZones();
 
         //- Check the surface zone definitions
-        void checkZones();
+        void checkZones(const bool verbose=true);
 
 
     // Modification
diff --git a/src/surfMesh/triSurface/triSurface.C b/src/surfMesh/triSurface/triSurface.C
index 51e2b52e1b2e595db53896b077b45f50e18f44b3..bf251788ee96353b53533f0217ee34bac8e8f421 100644
--- a/src/surfMesh/triSurface/triSurface.C
+++ b/src/surfMesh/triSurface/triSurface.C
@@ -637,16 +637,16 @@ void Foam::triSurface::swapPoints(pointField& pts)
 
 void Foam::triSurface::scalePoints(const scalar scaleFactor)
 {
-    // Avoid bad scaling
-    if (scaleFactor > SMALL && scaleFactor != 1.0)
+    // Avoid bad or no scaling
+    if (scaleFactor > SMALL && !equal(scaleFactor, 1))
     {
         // Remove all geometry dependent data
-        clearTopology();
+        this->clearTopology();
 
         // Adapt for new point positions
         MeshReference::movePoints(pointField());
 
-        storedPoints() *= scaleFactor;
+        this->storedPoints() *= scaleFactor;
     }
 }
 
@@ -666,6 +666,45 @@ void Foam::triSurface::cleanup(const bool verbose)
 }
 
 
+void Foam::triSurface::compactPoints(labelList& pointMap)
+{
+    this->clearOut();   // Topology changes
+
+    // Remove unused points while walking and renumbering faces
+    // in visit order - walk order as per localFaces()
+
+    labelList oldToCompact(this->points().size(), -1);
+    DynamicList<label> compactPointMap(oldToCompact.size());
+
+    for (auto& f : this->storedFaces())
+    {
+        for (label& pointi : f)
+        {
+            label compacti = oldToCompact[pointi];
+            if (compacti == -1)
+            {
+                compacti = compactPointMap.size();
+                oldToCompact[pointi] = compacti;
+                compactPointMap.append(pointi);
+            }
+            pointi = compacti;
+        }
+    }
+
+    pointField newPoints
+    (
+        UIndirectList<point>(this->points(), compactPointMap)
+    );
+
+    this->swapPoints(newPoints);
+
+    if (notNull(pointMap))
+    {
+        pointMap.transfer(compactPointMap);
+    }
+}
+
+
 Foam::List<Foam::surfZone>
 Foam::triSurface::sortedZones(labelList& faceMap) const
 {
diff --git a/src/surfMesh/triSurface/triSurface.H b/src/surfMesh/triSurface/triSurface.H
index 695f6fbe6818f0358fe021ff5b1e4e53858d67e8..770c08f8c4e2a6bcadb4950e7aa8e6f9a64ce280 100644
--- a/src/surfMesh/triSurface/triSurface.H
+++ b/src/surfMesh/triSurface/triSurface.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -397,18 +397,18 @@ public:
 
     // Access
 
-        const geometricSurfacePatchList& patches() const
+        const geometricSurfacePatchList& patches() const noexcept
         {
             return patches_;
         }
 
-        geometricSurfacePatchList& patches()
+        geometricSurfacePatchList& patches() noexcept
         {
             return patches_;
         }
 
         //- Return const access to the faces
-        inline const List<labelledTri>& surfFaces() const
+        const List<labelledTri>& surfFaces() const noexcept
         {
             return static_cast<const List<labelledTri>&>(*this);
         }
@@ -427,19 +427,19 @@ public:
 
 
         //- Face area vectors (normals)
-        inline const vectorField& Sf() const
+        const vectorField& Sf() const
         {
             return MeshReference::faceAreas();
         }
 
         //- Face area magnitudes
-        inline const scalarField& magSf() const
+        const scalarField& magSf() const
         {
             return MeshReference::magFaceAreas();
         }
 
         //- Face centres
-        inline const vectorField& Cf() const
+        const vectorField& Cf() const
         {
             return MeshReference::faceCentres();
         }
@@ -476,6 +476,14 @@ public:
         //- Remove non-valid triangles
         void cleanup(const bool verbose);
 
+        //- Remove unused points and renumber faces in local visit order
+        //
+        //  \param[out] pointMap from new to old points (optional)
+        void compactPoints
+        (
+            labelList& pointMap = const_cast<labelList&>(labelList::null())
+        );
+
         //- Fill faceZone with currentZone for every face reachable
         //  from facei without crossing edge marked in borderEdge.
         //  Note: faceZone has to be sized nFaces before calling this fun.