diff --git a/META-INFO/api-info b/META-INFO/api-info
index 1c8dcb8c5e52d5144f1fa851d5ae7cf9f6490b7e..e6c2adb53832eb9945aab63b3f49f29869acec0b 100644
--- a/META-INFO/api-info
+++ b/META-INFO/api-info
@@ -1,2 +1,2 @@
 api=2402
-patch=240220
+patch=240522
diff --git a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
index 7225bdcaaa8b23c8bf98c8ef4819ffcbeb4df207..3d5378c5774162ef4e40afb4f9d58db55efc2888 100644
--- a/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
+++ b/src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -310,7 +310,6 @@ void Foam::GAMGSolver::gatherMatrices
 
     const auto& boundaryMap = agglomeration_.boundaryMap(destLevel);
 
-    // Use PstreamBuffers
     PstreamBuffers pBufs(comm);
 
     // Send to master
@@ -327,35 +326,49 @@ void Foam::GAMGSolver::gatherMatrices
 
         const label proci = UPstream::myProcNo(comm);
 
-        labelList validInterface(interfaces.size(), -1);
+        // All interfaceBouCoeffs need to be sent across
+        bitSet validCoeffs(interfaces.size());
+        forAll(interfaceBouCoeffs, intI)
+        {
+            if (interfaceBouCoeffs.set(intI))
+            {
+                validCoeffs.set(intI);
+            }
+        }
+
+        // Only preserved interfaces need to be sent across
+        bitSet validInterface(interfaces.size());
         forAll(interfaces, intI)
         {
             const label allIntI = boundaryMap[proci][intI];
             if (interfaces.set(intI) && allIntI != -1)
             {
-                validInterface[intI] = intI;
+                validInterface.set(intI);
             }
         }
 
         UOPstream toMaster(UPstream::masterNo(), pBufs);
 
-        toMaster<< mat << token::SPACE << validInterface;
+        toMaster
+            << mat
+            << token::SPACE << validCoeffs
+            << token::SPACE << validInterface;
 
-        forAll(validInterface, intI)
+        for (const label intI : validCoeffs)
         {
-            if (validInterface[intI] != -1)
-            {
-                const auto& interface = refCast<const GAMGInterfaceField>
-                (
-                    interfaces[intI]
-                );
+            toMaster
+                << interfaceBouCoeffs[intI]
+                << interfaceIntCoeffs[intI];
+        }
+        for (const label intI : validInterface)
+        {
+            const auto& interface = refCast<const GAMGInterfaceField>
+            (
+                interfaces[intI]
+            );
 
-                toMaster
-                    << interfaceBouCoeffs[intI]
-                    << interfaceIntCoeffs[intI]
-                    << interface.type();
-                interface.write(toMaster);
-            }
+            toMaster << interface.type();
+            interface.write(toMaster);
         }
     }
 
@@ -371,10 +384,10 @@ void Foam::GAMGSolver::gatherMatrices
         lduInterfacePtrsList destInterfaces = destMesh.interfaces();
 
         // Master.
-        otherMats.setSize(nProcs-1);
-        otherBouCoeffs.setSize(nProcs-1);
-        otherIntCoeffs.setSize(nProcs-1);
-        otherInterfaces.setSize(nProcs-1);
+        otherMats.resize(nProcs-1);
+        otherBouCoeffs.resize(nProcs-1);
+        otherIntCoeffs.resize(nProcs-1);
+        otherInterfaces.resize(nProcs-1);
 
         for (const int proci : UPstream::subProcs(comm))
         {
@@ -384,60 +397,41 @@ void Foam::GAMGSolver::gatherMatrices
 
             otherMats.set(otherI, new lduMatrix(destMesh, fromProc));
 
-            // Receive number of/valid interfaces
-            // >= 0 : remote interface index
-            // -1   : invalid interface
-            const labelList validInterface(fromProc);
+            // Receive bit-sets of valid interfaceCoeffs/interfaces
+            const bitSet validCoeffs(fromProc);
+            const bitSet validInterface(fromProc);
 
-            otherBouCoeffs.set
-            (
-                otherI,
-                new FieldField<Field, scalar>(validInterface.size())
-            );
-            otherIntCoeffs.set
-            (
-                otherI,
-                new FieldField<Field, scalar>(validInterface.size())
-            );
-            otherInterfaces.set
-            (
-                otherI,
-                new PtrList<lduInterfaceField>(validInterface.size())
-            );
+            otherBouCoeffs.emplace_set(otherI, validCoeffs.size());
+            otherIntCoeffs.emplace_set(otherI, validCoeffs.size());
+            otherInterfaces.emplace_set(otherI, validInterface.size());
 
-            forAll(validInterface, intI)
+            // Receive individual interface contributions
+            for (const label intI : validCoeffs)
             {
-                if (validInterface[intI] != -1)
-                {
-                    otherBouCoeffs[otherI].set
-                    (
-                        intI,
-                        new scalarField(fromProc)
-                    );
-                    otherIntCoeffs[otherI].set
-                    (
-                        intI,
-                        new scalarField(fromProc)
-                    );
+                otherBouCoeffs[otherI].emplace_set(intI, fromProc);
+                otherIntCoeffs[otherI].emplace_set(intI, fromProc);
+            }
 
-                    const word coupleType(fromProc);
+            // Receive individual interface contributions
+            for (const label intI : validInterface)
+            {
+                const word coupleType(fromProc);
 
-                    const label allIntI = boundaryMap[proci][intI];
+                const label allIntI = boundaryMap[proci][intI];
 
-                    otherInterfaces[otherI].set
+                otherInterfaces[otherI].set
+                (
+                    intI,
+                    GAMGInterfaceField::New
                     (
-                        intI,
-                        GAMGInterfaceField::New
+                        coupleType,
+                        refCast<const GAMGInterface>
                         (
-                            coupleType,
-                            refCast<const GAMGInterface>
-                            (
-                                destInterfaces[allIntI]
-                            ),
-                            fromProc
-                        ).release()
-                    );
-                }
+                            destInterfaces[allIntI]
+                        ),
+                        fromProc
+                    ).release()
+                );
             }
         }
     }
diff --git a/src/fileFormats/stl/STLCore.C b/src/fileFormats/stl/STLCore.C
index 7e0c2d0ee1261985b6430935c888d71629d9d8d1..b77d68a7ad300d04583b9cb2571c8758ccb55454 100644
--- a/src/fileFormats/stl/STLCore.C
+++ b/src/fileFormats/stl/STLCore.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2023 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,6 +61,36 @@ static bool startsWithSolid(const char header[STLHeaderSize])
 }
 
 
+// Check if file size appears to be reasonable for an STL binary file.
+// Compare file size with that expected from number of tris
+// If this is not sensible, it may be an ASCII file
+//
+// sizeof(STLtriangle) = 50 bytes [int16 + 4 * (3 float)]
+
+inline static bool checkBinaryFileSize
+(
+    const int64_t nTris,
+    const Foam::fileName& file
+)
+{
+    // When checking the content size, account for the header size (80),
+    // but ignore the nTris information (int32_t) to give some rounding
+
+    const int64_t contentSize =
+    (
+        int64_t(Foam::fileSize(file))
+      - int64_t(STLHeaderSize)
+    );
+
+    return
+    (
+        (contentSize >= 0)
+     && (nTris >= contentSize/50)
+     && (nTris <= contentSize/25)
+    );
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 bool Foam::fileFormats::STLCore::isBinaryName
@@ -117,32 +147,31 @@ int Foam::fileFormats::STLCore::detectBinaryHeader
 
 
     // Read the number of triangles in the STL file
-    // (note: read as signed so we can check whether >2^31)
+    // (note: read as signed int so we can check whether >2^31).
+    //
+    // With nTris == 2^31, file size is 107.37 GB !
+    //
+    // However, the limit is more likely caused by the number of points
+    // that can be stored (label-size=32) when flattened for merging.
+    // So more like 715.8M triangles (~35.8 GB)
+
     int32_t nTris;
     is.read(reinterpret_cast<char*>(&nTris), sizeof(int32_t));
 
-    // Check that stream is OK and number of triangles is positive,
-    // if not this may be an ASCII file
-
-    bool bad = (!is || nTris < 0);
+    bool ok = (is && nTris >= 0);
 
-    if (!bad && unCompressed)
+    if (ok && unCompressed)
     {
-        // Compare file size with that expected from number of tris
-        // If this is not sensible, it may be an ASCII file
-
-        const off_t dataFileSize = Foam::fileSize(filename);
-
-        bad =
-            (
-                dataFileSize < STLHeaderSize
-             || nTris < (dataFileSize - STLHeaderSize)/50
-             || nTris > (dataFileSize - STLHeaderSize)/25
-            );
+        ok = checkBinaryFileSize(nTris, filename);
     }
 
+    //if (ok)
+    //{
+    //    InfoErr<< "stlb : " << nTris << " triangles" << nl;
+    //}
+
     // Return number of triangles if it appears to be BINARY and good.
-    return (bad ? 0 : nTris);
+    return (ok ? nTris : 0);
 }
 
 
@@ -190,32 +219,27 @@ Foam::fileFormats::STLCore::readBinaryHeader
             << exit(FatalError);
     }
 
-    // Read the number of triangles in the STl file
-    // (note: read as int so we can check whether >2^31)
+
+    // Read the number of triangles in the STL file
+    // (note: read as signed int so we can check whether >2^31).
+    //
+    // With nTris == 2^31, file size is 107.37 GB !
+    //
+    // However, the limit is more likely caused by the number of points
+    // that can be stored (label-size=32) when flattened for merging.
+    // So more like 715.8M triangles (~35.8 GB)
+
     int32_t nTris;
     is.read(reinterpret_cast<char*>(&nTris), sizeof(int32_t));
 
-    // Check that stream is OK and number of triangles is positive,
-    // if not this maybe an ASCII file
+    bool ok = (is && nTris >= 0);
 
-    bool bad = (!is || nTris < 0);
-
-    if (!bad && unCompressed)
+    if (ok && unCompressed)
     {
-        // Compare file size with that expected from number of tris
-        // If this is not sensible, it may be an ASCII file
-
-        const off_t dataFileSize = Foam::fileSize(filename);
-
-        bad =
-            (
-                dataFileSize < STLHeaderSize
-             || nTris < (dataFileSize - STLHeaderSize)/50
-             || nTris > (dataFileSize - STLHeaderSize)/25
-            );
+        ok = checkBinaryFileSize(nTris, filename);
     }
 
-    if (bad)
+    if (!ok)
     {
         FatalErrorInFunction
             << "problem reading number of triangles, perhaps file is not binary"
diff --git a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
index 6fa28be7ea89a7cc15d4918b963b99d3350549bd..3e6c9b85c344aae4f4e94e6d53754f07ef3eb125 100644
--- a/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/sampledMeshedSurface/sampledMeshedSurfaceTemplates.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,7 +39,11 @@ Foam::sampledMeshedSurface::sampleOnFaces
 {
     const Type deflt
     (
-        defaultValues_.getOrDefault<Type>(sampler.psi().name(), Zero)
+        defaultValues_.getOrDefault<Type>
+        (
+            sampler.psi().name(),
+            Foam::zero{}
+        )
     );
 
     const labelList& elements = sampleElements_;
@@ -71,13 +75,16 @@ Foam::sampledMeshedSurface::sampleOnFaces
 
     const polyBoundaryMesh& pbm = mesh().boundaryMesh();
 
-    Field<Type> bVals(mesh().nBoundaryFaces(), Zero);
+    Field<Type> bVals(mesh().nBoundaryFaces(), deflt);
 
     const auto& bField = sampler.psi().boundaryField();
 
     forAll(bField, patchi)
     {
-        SubList<Type>(bVals, pbm[patchi].range()) = bField[patchi];
+        // Note: restrict transcribing to actual size of the patch field
+        // - handles "empty" patch type etc.
+        const auto& pfld = bField[patchi];
+        SubList<Type>(bVals, pfld.size(), pbm[patchi].offset()) = pfld;
     }
 
     // Sample within the flat boundary field
@@ -109,7 +116,11 @@ Foam::sampledMeshedSurface::sampleOnPoints
 {
     const Type deflt
     (
-        defaultValues_.getOrDefault<Type>(interpolator.psi().name(), Zero)
+        defaultValues_.getOrDefault<Type>
+        (
+            interpolator.psi().name(),
+            Foam::zero{}
+        )
     );
 
     const labelList& elements = sampleElements_;