diff --git a/applications/test/globalMeshData/Make/files b/applications/test/globalMeshData/Make/files
index 80e30e0d4a5e9d3d0fc57b598badcad46e36fc3e..cb776b3445c2eb0f1abc9a1f78056eb55bcbe43b 100644
--- a/applications/test/globalMeshData/Make/files
+++ b/applications/test/globalMeshData/Make/files
@@ -1,3 +1,3 @@
-Test-globalMeshData.C
+Test-globalMeshData.cxx
 
 EXE = $(FOAM_USER_APPBIN)/Test-globalMeshData
diff --git a/applications/test/globalMeshData/Test-globalMeshData.C b/applications/test/globalMeshData/Test-globalMeshData.cxx
similarity index 97%
rename from applications/test/globalMeshData/Test-globalMeshData.C
rename to applications/test/globalMeshData/Test-globalMeshData.cxx
index ffc357a5493c74eca0d20014439ea6c0c8b04564..b078620065453714333b4690d3ede9e8d5cb9cc8 100644
--- a/applications/test/globalMeshData/Test-globalMeshData.C
+++ b/applications/test/globalMeshData/Test-globalMeshData.cxx
@@ -182,18 +182,11 @@ int main(int argc, char *argv[])
         const labelListList& transformedSlaves =
             globalData.globalPointTransformedBoundaryFaces();
 
-        const label nBnd = mesh.nBoundaryFaces();
-
         pointField fc(globalPointBoundaryFacesMap.constructSize());
-        SubList<point>(fc, nBnd) =
+        pointField::subList(fc, mesh.nBoundaryFaces()) =
             primitivePatch
             (
-                SubList<face>
-                (
-                    mesh.faces(),
-                    nBnd,
-                    mesh.nInternalFaces()
-                ),
+                mesh.boundaryMesh().faces(),
                 mesh.points()
             ).faceCentres();
 
diff --git a/applications/test/syncTools/Make/files b/applications/test/syncTools/Make/files
index a9c9a0e01367398b6e11ceb0b50793838fd5b485..a1ddebcb63db624a0efaaac4814b3c5830570c96 100644
--- a/applications/test/syncTools/Make/files
+++ b/applications/test/syncTools/Make/files
@@ -1,3 +1,3 @@
-Test-syncTools.C
+Test-syncTools.cxx
 
 EXE = $(FOAM_USER_APPBIN)/Test-syncTools
diff --git a/applications/test/syncTools/Test-syncTools.C b/applications/test/syncTools/Test-syncTools.cxx
similarity index 99%
rename from applications/test/syncTools/Test-syncTools.C
rename to applications/test/syncTools/Test-syncTools.cxx
index 70c288274d5bfdf533ffa36f8dffa1eafb8e8575..05bcbf46e269ce957222d4cacdc7be8a259aad95 100644
--- a/applications/test/syncTools/Test-syncTools.C
+++ b/applications/test/syncTools/Test-syncTools.cxx
@@ -195,16 +195,7 @@ void testSparseData(const polyMesh& mesh, Random& rndGen)
         << "Position test of sparse data only correct for cases without cyclics"
         << " with shared points." << endl;
 
-    primitivePatch allBoundary
-    (
-        SubList<face>
-        (
-            mesh.faces(),
-            mesh.nBoundaryFaces(),
-            mesh.nInternalFaces()
-        ),
-        mesh.points()
-    );
+    primitivePatch allBoundary(mesh.boundaryMesh().faces(), mesh.points());
     const pointField& localPoints = allBoundary.localPoints();
 
 
diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
index 34f157dfa0ad5b922b631c043ee5088f7fee6c0e..cdce544196eae85d057f6f2d1c49d88b92fa10ff 100644
--- a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
+++ b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
@@ -417,14 +417,7 @@ int main(int argc, char *argv[])
 
 
     // Get calculating engine for all of outside
-    const SubList<face> outsideFaces
-    (
-        mesh.faces(),
-        mesh.nBoundaryFaces(),
-        mesh.nInternalFaces()
-    );
-
-    primitivePatch allBoundary(outsideFaces, mesh.points());
+    primitivePatch allBoundary(mesh.boundaryMesh().faces(), mesh.points());
 
 
     // Look up mesh labels and convert to input for boundaryCutter.
diff --git a/applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H b/applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H
index e46205b7bcceffdab20ac1ec23f396a8e66d56ea..b35c2a76c976931aaf47a50b7450cf558781a032 100644
--- a/applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H
+++ b/applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H
@@ -104,18 +104,13 @@ const polyMesh& topoMesh = topoMeshPtr();
 
 // Block boundary faces
 {
-    const label nIntFaces = topoMesh.nInternalFaces();
+    const auto& patches = topoMesh.boundaryMesh();
     const label nBndFaces = topoMesh.nBoundaryFaces();
 
-    faceList bndFaces
-    (
-        faceList::subList(topoMesh.faces(), nBndFaces, nIntFaces)
-    );
-
     vtk::surfaceWriter writer
     (
         topoMesh.points(),
-        bndFaces,
+        faceList(patches.faces()),
         vtk::formatType::INLINE_ASCII,
         runTime.path()/"blockFaces"
     );
@@ -132,7 +127,7 @@ const polyMesh& topoMesh = topoMeshPtr();
 
         for (const polyPatch& pp : patches)
         {
-            label bndFacei = pp.start() - nIntFaces;
+            label bndFacei = pp.offset();
             label meshFacei = pp.start();
 
             forAll(pp, bfacei)
diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index 4a93e7117b3040f591cbaa6bd19c20f7ad695ec3..364613e238faff4f09a52f5cc5a7970def7a9d07 100644
--- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -608,12 +608,7 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
 {
     primitivePatch tmpBoundaryFaces
     (
-        SubList<face>
-        (
-            mesh_.faces(),
-            mesh_.nBoundaryFaces(),
-            mesh_.nInternalFaces()
-        ),
+        mesh_.boundaryMesh().faces(),
         mesh_.points()
     );
 
diff --git a/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C b/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C
index 971e5fa169579a7bec4d027f3b78d427068c76d2..14eebe1ec843e39a44c9bc47045bb378b92e90f4 100644
--- a/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C
+++ b/applications/utilities/mesh/manipulation/polyDualMesh/polyDualMeshApp.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016,2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -137,16 +137,7 @@ void simpleMarkFeatures
     // multiple 'half' cells.
 
     // Addressing for all outside faces
-    primitivePatch allBoundary
-    (
-        SubList<face>
-        (
-            mesh.faces(),
-            mesh.nBoundaryFaces(),
-            mesh.nInternalFaces()
-        ),
-        mesh.points()
-    );
+    primitivePatch allBoundary(patches.faces(), mesh.points());
 
     // Check for non-manifold points (surface pinched at point)
     allBoundary.checkPointManifold(false, &singleCellFeaturePointSet);
@@ -240,7 +231,9 @@ void simpleMarkFeatures
     // ~~~~~~~~~~~~~~~~~~~~~~~~~
 
     // Face centres that need inclusion in the dual mesh
-    labelHashSet featureFaceSet(mesh.nBoundaryFaces());
+    labelHashSet featureFaceSet;
+    featureFaceSet.reserve(mesh.nBoundaryFaces());
+
     // A. boundary faces.
     for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++)
     {
diff --git a/src/OpenFOAM/meshes/pointMesh/pointPatches/meshPointPatch/meshPointPatch.C b/src/OpenFOAM/meshes/pointMesh/pointPatches/meshPointPatch/meshPointPatch.C
index 2e9c845c75cd2eacec46fd6aa351bed7b4938b4f..c86861ac9272a07c80da3afeeb70d81f7fa3d552 100644
--- a/src/OpenFOAM/meshes/pointMesh/pointPatches/meshPointPatch/meshPointPatch.C
+++ b/src/OpenFOAM/meshes/pointMesh/pointPatches/meshPointPatch/meshPointPatch.C
@@ -160,14 +160,7 @@ void Foam::meshPointPatch::movePoints(PstreamBuffers&, const pointField& p)
     // Recalculate the point normals? Something like
     //if (owner())
     //{
-    //    const SubList<face> subFaces
-    //    (
-    //        mesh.faces(),
-    //        mesh.nBoundaryFaces(),
-    //        mesh.nInternalFaces()
-    //    );
-    //    const primitivePatch pp(subFaces, mesh.points());
-    //
+    //    const primitivePatch pp(mesh.boundaryMesh().faces(), mesh.points());
     //
     //    for (const label pointi : meshPoints())
     //    {
@@ -180,9 +173,8 @@ void Foam::meshPointPatch::movePoints(PstreamBuffers&, const pointField& p)
     //        const auto& point
     //
     //
-
 }
-    
+
 
 void Foam::meshPointPatch::updateMesh(PstreamBuffers&)
 {
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
index db304f4b2ec840287e5ec7888e7c26c77fd42cd9..ca3c41504ead451c92a858464a5e2190600d8f33 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2024 OpenCFD Ltd.
+    Copyright (C) 2018-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -91,7 +91,7 @@ Foam::polyPatch::polyPatch
     patchIdentifier(name, index),
     primitivePatch
     (
-        faceSubList(bm.mesh().faces(), size, start),
+        SubList<face>(bm.mesh().faces(), size, start),
         bm.mesh().points()
     ),
     start_(start),
@@ -118,7 +118,7 @@ Foam::polyPatch::polyPatch
     patchIdentifier(name, index, physicalType, inGroups),
     primitivePatch
     (
-        faceSubList(bm.mesh().faces(), size, start),
+        SubList<face>(bm.mesh().faces(), size, start),
         bm.mesh().points()
     ),
     start_(start),
@@ -138,7 +138,7 @@ Foam::polyPatch::polyPatch
     patchIdentifier(name, dict, index),
     primitivePatch
     (
-        faceSubList
+        SubList<face>
         (
             bm.mesh().faces(),
             dict.get<label>("nFaces"),
@@ -165,7 +165,7 @@ Foam::polyPatch::polyPatch
     patchIdentifier(pp),
     primitivePatch
     (
-        faceSubList
+        SubList<face>
         (
             bm.mesh().faces(),
             pp.size(),
@@ -190,7 +190,7 @@ Foam::polyPatch::polyPatch
     patchIdentifier(pp, index),
     primitivePatch
     (
-        faceSubList
+        SubList<face>
         (
             bm.mesh().faces(),
             newSize,
@@ -215,7 +215,7 @@ Foam::polyPatch::polyPatch
     patchIdentifier(pp, index),
     primitivePatch
     (
-        faceSubList
+        SubList<face>
         (
             bm.mesh().faces(),
             mapAddressing.size(),
@@ -305,6 +305,25 @@ const Foam::polyBoundaryMesh& Foam::polyPatch::boundaryMesh() const noexcept
 }
 
 
+const Foam::faceList::subList Foam::polyPatch::faces() const
+{
+    return patchSlice(boundaryMesh().mesh().faces());
+}
+
+
+const Foam::labelList::subList Foam::polyPatch::faceOwner() const
+{
+    return patchSlice(boundaryMesh().mesh().faceOwner());
+}
+
+
+// Potentially useful to simplify logic elsewhere?
+// const Foam::labelList::subList Foam::polyPatch::faceNeighbour() const
+// {
+//     return labelList::subList();
+// }
+
+
 const Foam::vectorField::subField Foam::polyPatch::faceCentres() const
 {
     return patchSlice(boundaryMesh().mesh().faceCentres());
diff --git a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
index 02f9901e20fc8210f00ed43962ab3a285a9a0939..2b6d20eccd65c0aa047a47a8d2c435075dbb8bbd 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
+++ b/src/OpenFOAM/meshes/polyMesh/polyPatches/polyPatch/polyPatch.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2024 OpenCFD Ltd.
+    Copyright (C) 2015-2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -435,6 +435,12 @@ public:
 
         // Geometric data; point list required
 
+            //- Return mesh faces for the patch
+            const faceList::subList faces() const;
+
+            //- Return face owner for the patch
+            const labelList::subList faceOwner() const;
+
             //- Return face centres
             const vectorField::subField faceCentres() const;
 
diff --git a/src/conversion/polyDualMesh/polyDualMesh.C b/src/conversion/polyDualMesh/polyDualMesh.C
index f60f0106d7c42a99b4b4f27333dcc534a0035a87..72807e8f16bba8e4c7cdfb045669c9534102b4ca 100644
--- a/src/conversion/polyDualMesh/polyDualMesh.C
+++ b/src/conversion/polyDualMesh/polyDualMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
+    Copyright (C) 2022,2025 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -723,16 +723,7 @@ void Foam::polyDualMesh::calcDual
 
 
     // Get patch for all of outside
-    primitivePatch allBoundary
-    (
-        SubList<face>
-        (
-            mesh.faces(),
-            mesh.nBoundaryFaces(),
-            mesh.nInternalFaces()
-        ),
-        mesh.points()
-    );
+    primitivePatch allBoundary(patches.faces(), mesh.points());
 
     // Correspondence from patch edge to mesh edge.
     labelList meshEdges
@@ -1467,31 +1458,19 @@ void Foam::polyDualMesh::calcFeatures
     labelList& featurePoints
 )
 {
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
     // Create big primitivePatch for all outside.
-    primitivePatch allBoundary
-    (
-        SubList<face>
-        (
-            mesh.faces(),
-            mesh.nBoundaryFaces(),
-            mesh.nInternalFaces()
-        ),
-        mesh.points()
-    );
+    primitivePatch allBoundary(patches.faces(), mesh.points());
 
     // For ease of use store patch number per face in allBoundary.
     labelList allRegion(allBoundary.size());
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
     forAll(patches, patchi)
     {
         const polyPatch& pp = patches[patchi];
 
-        forAll(pp, i)
-        {
-            allRegion[i + pp.start() - mesh.nInternalFaces()] = patchi;
-        }
+        allRegion.slice(pp.offset(), pp.size()) = patchi;
     }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/IntegralScaleBox/IntegralScaleBox.C b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/IntegralScaleBox/IntegralScaleBox.C
index 47cebe3c6a8747cfdc421d5dd98f280d5cc63b24..3e0077a2093f27511eb232e015111489c3b9e7c3 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/IntegralScaleBox/IntegralScaleBox.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/turbulentDigitalFilterInlet/IntegralScaleBox/IntegralScaleBox.C
@@ -372,7 +372,7 @@ void Foam::turbulence::IntegralScaleBox<Type>::calcPatch()
         (
             new primitivePatch
             (
-                SubList<face>(patchFaces_, patchFaces_.size()),
+                SubList<face>(patchFaces_),
                 patchPoints_
             )
         );
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C
index 31f4f9f0a47f839a3111953a55ec5b96d93221e4..699a85116e343459e5bd416a36bbadb9ef162f72 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolation.C
@@ -71,22 +71,12 @@ void Foam::volPointInterpolation::calcBoundaryAddressing()
             << endl;
     }
 
-    boundaryPtr_.reset
-    (
-        new primitivePatch
-        (
-            SubList<face>
-            (
-                mesh().faces(),
-                mesh().nBoundaryFaces(),
-                mesh().nInternalFaces()
-            ),
-            mesh().points()
-        )
-    );
-    const primitivePatch& boundary = boundaryPtr_();
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+    boundaryPtr_.reset(new primitivePatch(pbm.faces(), mesh().points()));
+    const auto& boundary = *boundaryPtr_;
 
-    boundaryIsPatchFace_.setSize(boundary.size());
+    boundaryIsPatchFace_.resize_nocopy(boundary.size());
     boundaryIsPatchFace_ = false;
 
     // Store per mesh point whether it is on any 'real' patch. Currently
@@ -94,8 +84,6 @@ void Foam::volPointInterpolation::calcBoundaryAddressing()
     // bitSet. Tbd)
     boolList isPatchPoint(mesh().nPoints(), false);
 
-    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-
     // Get precalculated volField only so we can use coupled() tests for
     // cyclicAMI
     const surfaceScalarField& magSf = mesh().magSf();
@@ -110,18 +98,11 @@ void Foam::volPointInterpolation::calcBoundaryAddressing()
          && !magSf.boundaryField()[patchi].coupled()
         )
         {
-            label bFacei = pp.start()-mesh().nInternalFaces();
+            boundaryIsPatchFace_.set(labelRange(pp.offset(), pp.size()));
 
-            forAll(pp, i)
+            for (const auto& f : pp.faces())
             {
-                boundaryIsPatchFace_[bFacei] = true;
-
-                const face& f = boundary[bFacei++];
-
-                forAll(f, fp)
-                {
-                    isPatchPoint[f[fp]] = true;
-                }
+                UIndirectList<bool>(isPatchPoint, f) = true;
             }
         }
     }
@@ -136,27 +117,15 @@ void Foam::volPointInterpolation::calcBoundaryAddressing()
     );
 
     // Convert to bitSet
-    isPatchPoint_.setSize(mesh().nPoints());
+    isPatchPoint_.reset();
+    isPatchPoint_.resize(mesh().nPoints());
     isPatchPoint_.assign(isPatchPoint);
 
     if (debug)
     {
-        label nPatchFace = 0;
-        forAll(boundaryIsPatchFace_, i)
-        {
-            if (boundaryIsPatchFace_[i])
-            {
-                nPatchFace++;
-            }
-        }
-        label nPatchPoint = 0;
-        forAll(isPatchPoint_, i)
-        {
-            if (isPatchPoint_[i])
-            {
-                nPatchPoint++;
-            }
-        }
+        label nPatchFace = boundaryIsPatchFace_.count();
+        label nPatchPoint = isPatchPoint_.count();
+
         Pout<< "boundary:" << nl
             << "    faces :" << boundary.size() << nl
             << "    of which on proper patch:" << nPatchFace << nl
diff --git a/src/lumpedPointMotion/movement/lumpedPointMovementWriter.C b/src/lumpedPointMotion/movement/lumpedPointMovementWriter.C
index 84cd47ada606942f652f8b76257a56c17e8e578b..46d293cb2a32831ff5f6d2e5b2eccbc6f284e08d 100644
--- a/src/lumpedPointMotion/movement/lumpedPointMovementWriter.C
+++ b/src/lumpedPointMotion/movement/lumpedPointMovementWriter.C
@@ -235,11 +235,7 @@ void Foam::lumpedPointMovement::writeZonesVTP
     {
         const labelList& faceToPoint = patchControls_[patchi].faceToPoint_;
 
-        primitivePatch pp
-        (
-            SubList<face>(mesh.faces(), patches[patchi].range()),
-            points0
-        );
+        primitivePatch pp(patches[patchi].faces(), points0);
 
         writer.piece(pp.localPoints(), pp.localFaces());
 
diff --git a/src/optimisation/adjointOptimisation/adjoint/interpolation/volPointInterpolation/volPointInterpolationAdjoint.C b/src/optimisation/adjointOptimisation/adjoint/interpolation/volPointInterpolation/volPointInterpolationAdjoint.C
index 0d63d8cff345d7562d95a7bed048c6cda7c50b93..a9156ba1fda54a43cc8cada82909d1fb8c7b9bd3 100644
--- a/src/optimisation/adjointOptimisation/adjoint/interpolation/volPointInterpolation/volPointInterpolationAdjoint.C
+++ b/src/optimisation/adjointOptimisation/adjoint/interpolation/volPointInterpolation/volPointInterpolationAdjoint.C
@@ -55,22 +55,12 @@ void Foam::volPointInterpolationAdjoint::calcBoundaryAddressing()
             << endl;
     }
 
-    boundaryPtr_.reset
-    (
-        new primitivePatch
-        (
-            SubList<face>
-            (
-                mesh().faces(),
-                mesh().nBoundaryFaces(),
-                mesh().nInternalFaces()
-            ),
-            mesh().points()
-        )
-    );
-    const primitivePatch& boundary = boundaryPtr_();
+    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+    boundaryPtr_.reset(new primitivePatch(pbm.faces(), mesh().points()));
+    const auto& boundary = *boundaryPtr_;
 
-    boundaryIsPatchFace_.setSize(boundary.size());
+    boundaryIsPatchFace_.resize_nocopy(boundary.size());
     boundaryIsPatchFace_ = false;
 
     // Store per mesh point whether it is on any 'real' patch. Currently
@@ -79,8 +69,6 @@ void Foam::volPointInterpolationAdjoint::calcBoundaryAddressing()
     boolList isPatchPoint(mesh().nPoints(), false);
     boolList isSymmetryPoint(mesh().nPoints(), false);
 
-    const polyBoundaryMesh& pbm = mesh().boundaryMesh();
-
     // Get precalculated volField only so we can use coupled() tests for
     // cyclicAMI
     const surfaceScalarField& magSf = mesh().magSf();
@@ -97,27 +85,16 @@ void Foam::volPointInterpolationAdjoint::calcBoundaryAddressing()
          && !isA<symmetryPlanePolyPatch>(pp)
         )
         {
-            label bFacei = pp.start()-mesh().nInternalFaces();
+            boundaryIsPatchFace_.set(labelRange(pp.offset(), pp.size()));
 
-            forAll(pp, i)
+            for (const auto& f : pp.faces())
             {
-                boundaryIsPatchFace_[bFacei] = true;
-
-                const face& f = boundary[bFacei++];
-
-                forAll(f, fp)
-                {
-                    isPatchPoint[f[fp]] = true;
-                }
+                UIndirectList<bool>(isPatchPoint, f) = true;
             }
         }
         else if (isA<symmetryPolyPatch>(pp) || isA<symmetryPlanePolyPatch>(pp))
         {
-            const labelList& meshPoints = pp.meshPoints();
-            for (const label pointI : meshPoints)
-            {
-                isSymmetryPoint[pointI] = true;
-            }
+            UIndirectList<bool>(isSymmetryPoint, pp.meshPoints()) = true;
         }
     }
 
@@ -131,30 +108,19 @@ void Foam::volPointInterpolationAdjoint::calcBoundaryAddressing()
     );
 
     // Convert to bitSet
-    isPatchPoint_.setSize(mesh().nPoints());
+    isPatchPoint_.reset();
+    isPatchPoint_.resize(mesh().nPoints());
     isPatchPoint_.assign(isPatchPoint);
 
-    isSymmetryPoint_.setSize(mesh().nPoints());
+    isSymmetryPoint_.reset();
+    isSymmetryPoint_.resize(mesh().nPoints());
     isSymmetryPoint_.assign(isSymmetryPoint);
 
     if (debug)
     {
-        label nPatchFace = 0;
-        forAll(boundaryIsPatchFace_, i)
-        {
-            if (boundaryIsPatchFace_[i])
-            {
-                nPatchFace++;
-            }
-        }
-        label nPatchPoint = 0;
-        forAll(isPatchPoint_, i)
-        {
-            if (isPatchPoint_[i])
-            {
-                nPatchPoint++;
-            }
-        }
+        label nPatchFace = boundaryIsPatchFace_.count();
+        label nPatchPoint = isPatchPoint_.count();
+
         Pout<< "boundary:" << nl
             << "    faces :" << boundary.size() << nl
             << "    of which on proper patch:" << nPatchFace << nl
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index 30516a5518ca18721cf5dd0f09b544e90a88cebe..64b78b04c7b1379ce219db4ddc5123aa34e0c1e7 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.C
@@ -375,46 +375,20 @@ Foam::MeshedSurface<Face>::MeshedSurface(const surfMesh& mesh)
 template<class Face>
 Foam::MeshedSurface<Face>::MeshedSurface
 (
-    const polyBoundaryMesh& bMesh,
+    const polyBoundaryMesh& pbm,
     const bool useGlobalPoints
 )
 :
     MeshedSurface<Face>()
 {
-    const polyMesh& mesh = bMesh.mesh();
-    const polyPatchList& bPatches = bMesh;
-
-    // Get a single patch for all boundaries
-    primitivePatch allBoundary
-    (
-        SubList<face>
-        (
-            mesh.faces(),
-            mesh.nBoundaryFaces(),
-            mesh.nInternalFaces()
-        ),
-        mesh.points()
-    );
-
-    // use global/local points:
-    const pointField& bPoints =
-    (
-        useGlobalPoints ? mesh.points() : allBoundary.localPoints()
-    );
-
-    // global/local face addressing:
-    const List<Face>& bFaces =
-    (
-        useGlobalPoints ? allBoundary : allBoundary.localFaces()
-    );
-
+    const polyMesh& mesh = pbm.mesh();
 
     // create zone list
-    surfZoneList newZones(bPatches.size());
+    surfZoneList newZones(pbm.size());
 
     label startFacei = 0;
     label nZone = 0;
-    for (const polyPatch& p : bPatches)
+    for (const polyPatch& p : pbm)
     {
         if (p.size())
         {
@@ -433,8 +407,27 @@ Foam::MeshedSurface<Face>::MeshedSurface
 
     newZones.setSize(nZone);
 
+    // Get a single patch for all boundaries
+    primitivePatch allBoundary(pbm.faces(), mesh.points());
+
     // Face type as per polyBoundaryMesh
-    MeshedSurface<face> surf(bPoints, bFaces, newZones);
+    MeshedSurface<face> surf;
+
+    if (useGlobalPoints)
+    {
+        // use global points, global face addressing
+        surf = MeshedSurface<face>(mesh.points(), allBoundary, newZones);
+    }
+    else
+    {
+        // use local points, local face addressing
+        surf = MeshedSurface<face>
+        (
+            allBoundary.localPoints(),
+            allBoundary.localFaces(),
+            newZones
+        );
+    }
 
     this->transcribe(surf);
 }
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.H b/src/surfMesh/MeshedSurface/MeshedSurface.H
index 8f7214257e34a276e6a2cf419117d29dc01e85fa..81be8efbe3e5c9c36eda223e1eec9e71d2a95560 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.H
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.H
@@ -288,7 +288,7 @@ public:
         //- Construct from a boundary mesh with local points/faces
         MeshedSurface
         (
-            const polyBoundaryMesh& bMesh,
+            const polyBoundaryMesh& pbm,
             const bool globalPoints = false
         );