BUG: GAMG: processorAgglomeration not merging proc coeffs. Fixes #3172

- requires a field size check when copying into the flat boundary

META-INFO/api-info

@@ -1,2 +1,2 @@
 api=2312
-patch=0
+patch=240220

etc/config.csh/setup

@@ -197,8 +197,8 @@ if ( "$FOAM_MPI" != dummy ) then
     _foamAddLib "${FOAM_LIBBIN}/${FOAM_MPI}"
 endif
 
-# OpenFOAM user, group libraries
-_foamAddLib  "${FOAM_USER_LIBBIN}:${FOAM_SITE_LIBBIN}"
+_foamAddLib  "$FOAM_SITE_LIBBIN"  # OpenFOAM group libraries
+_foamAddLib  "$FOAM_USER_LIBBIN"  # OpenFOAM user  libraries
 
 if ( -d "$WM_PROJECT_DIR/doc/man1" ) then
     _foamAddMan "$WM_PROJECT_DIR/doc"

etc/config.sh/setup

@@ -229,8 +229,8 @@ then
     _foamAddLib "$FOAM_LIBBIN/$FOAM_MPI"
 fi
 
-# OpenFOAM user, group libraries
-_foamAddLib  "$FOAM_USER_LIBBIN:$FOAM_SITE_LIBBIN"
+_foamAddLib  "$FOAM_SITE_LIBBIN"  # OpenFOAM group libraries
+_foamAddLib  "$FOAM_USER_LIBBIN"  # OpenFOAM user  libraries
 
 if [ -d "$WM_PROJECT_DIR/doc/man1" ]
 then

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.
@@ -332,35 +332,47 @@ 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);
         }
     }
 
@@ -389,60 +401,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 bitSet 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()
+                );
             }
         }
     }

src/OpenFOAM/matrices/schemes/schemesLookup.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -133,7 +133,7 @@ Foam::schemesLookup::schemesLookup
             obr.time().system(),
             obr,
             rOpt,
-            IOobject::NO_WRITE
+            IOobjectOption::NO_WRITE
         ),
         fallback
     ),
@@ -153,18 +153,23 @@ Foam::schemesLookup::schemesLookup
     fluxRequiredDefault_(false),
     steady_(false)
 {
-    // Treat as MUST_READ_IF_MODIFIED whenever possible
+    // Treat as READ_MODIFIED whenever possible
     if
     (
-        readOpt() == IOobject::MUST_READ
+        readOpt() == IOobjectOption::MUST_READ
      || (isReadOptional() && headerOk())
     )
     {
-        readOpt(IOobject::MUST_READ_IF_MODIFIED);
+        readOpt(IOobjectOption::READ_MODIFIED);
         addWatch();
     }
 
-    if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+    // Update: from values read or copied in
+    if
+    (
+        readOpt() == IOobjectOption::READ_MODIFIED
+     || !dictionary::empty()
+    )
     {
         read(selectedDict());
     }

src/OpenFOAM/matrices/solution/solution.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -174,7 +174,7 @@ Foam::solution::solution
             obr.time().system(),
             obr,
             rOpt,
-            IOobject::NO_WRITE
+            IOobjectOption::NO_WRITE
         ),
         fallback
     ),
@@ -184,18 +184,23 @@ Foam::solution::solution
     eqnRelaxDict_(),
     solvers_()
 {
-    // Treat as MUST_READ_IF_MODIFIED whenever possible
+    // Treat as READ_MODIFIED whenever possible
     if
     (
-        readOpt() == IOobject::MUST_READ
+        readOpt() == IOobjectOption::MUST_READ
      || (isReadOptional() && headerOk())
     )
     {
-        readOpt(IOobject::MUST_READ_IF_MODIFIED);
+        readOpt(IOobjectOption::READ_MODIFIED);
         addWatch();
     }
 
-    if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+    // Update: from values read or copied in
+    if
+    (
+        readOpt() == IOobjectOption::READ_MODIFIED
+     || !dictionary::empty()
+    )
     {
         read(selectedDict());
     }

src/finiteVolume/fields/fvPatchFields/constraint/cyclicACMI/cyclicACMIFvPatchField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -444,7 +444,14 @@ void Foam::cyclicACMIFvPatchField<Type>::initEvaluate
             << " starting send&receive"
             << endl;
 
-        if (!this->ready())
+        // Bypass polyPatch to get nbrId.
+        // - use cyclicACMIFvPatch::neighbPatch() virtual instead
+        const cyclicACMIFvPatch& neighbPatch = cyclicACMIPatch_.neighbPatch();
+        const labelUList& nbrFaceCells = neighbPatch.faceCells();
+        const Field<Type> pnf(this->primitiveField(), nbrFaceCells);
+
+        // Assert that all receives are known to have finished
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -453,11 +460,8 @@ void Foam::cyclicACMIFvPatchField<Type>::initEvaluate
                 << abort(FatalError);
         }
 
-        // By-pass polyPatch to get nbrId. Instead use cyclicACMIFvPatch virtual
-        // neighbPatch()
-        const cyclicACMIFvPatch& neighbPatch = cyclicACMIPatch_.neighbPatch();
-        const labelUList& nbrFaceCells = neighbPatch.faceCells();
-        const Field<Type> pnf(this->primitiveField(), nbrFaceCells);
+        // Assume that sends are also OK
+        sendRequests_.clear();
 
         cyclicACMIPatch_.initInterpolate
         (
@@ -515,6 +519,10 @@ void Foam::cyclicACMIFvPatchField<Type>::evaluate
             ).ptr()
         );
 
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
+
+
         auto& patchNeighbourField = patchNeighbourFieldPtr_.ref();
 
         if (doTransform())
@@ -559,7 +567,16 @@ void Foam::cyclicACMIFvPatchField<Type>::initInterfaceMatrixUpdate
             << " starting send&receive"
             << endl;
 
-        if (!this->ready())
+        const labelUList& nbrFaceCells =
+            lduAddr.patchAddr(cyclicACMIPatch_.neighbPatchID());
+
+        solveScalarField pnf(psiInternal, nbrFaceCells);
+
+        // Transform according to the transformation tensors
+        transformCoupleField(pnf, cmpt);
+
+        // Assert that all receives are known to have finished
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -568,13 +585,8 @@ void Foam::cyclicACMIFvPatchField<Type>::initInterfaceMatrixUpdate
                 << abort(FatalError);
         }
 
-        const labelUList& nbrFaceCells =
-            lduAddr.patchAddr(cyclicACMIPatch_.neighbPatchID());
-
-        solveScalarField pnf(psiInternal, nbrFaceCells);
-
-        // Transform according to the transformation tensors
-        transformCoupleField(pnf, cmpt);
+        // Assume that sends are also OK
+        sendRequests_.clear();
 
         cyclicACMIPatch_.initInterpolate
         (
@@ -635,6 +647,9 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
                 recvRequests_,
                 scalarRecvBufs_
             );
+
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
     }
     else
     {
@@ -676,7 +691,16 @@ void Foam::cyclicACMIFvPatchField<Type>::initInterfaceMatrixUpdate
                 << exit(FatalError);
         }
 
-        if (!this->ready())
+        const labelUList& nbrFaceCells =
+            lduAddr.patchAddr(cyclicACMIPatch_.neighbPatchID());
+
+        Field<Type> pnf(psiInternal, nbrFaceCells);
+
+        // Transform according to the transformation tensors
+        transformCoupleField(pnf);
+
+        // Assert that all receives are known to have finished
+        if (!recvRequests_.empty())
         {
             FatalErrorInFunction
                 << "Outstanding recv request(s) on patch "
@@ -685,13 +709,8 @@ void Foam::cyclicACMIFvPatchField<Type>::initInterfaceMatrixUpdate
                 << abort(FatalError);
         }
 
-        const labelUList& nbrFaceCells =
-            lduAddr.patchAddr(cyclicACMIPatch_.neighbPatchID());
-
-        Field<Type> pnf(psiInternal, nbrFaceCells);
-
-        // Transform according to the transformation tensors
-        transformCoupleField(pnf);
+        // Assume that sends are also OK
+        sendRequests_.clear();
 
         cyclicACMIPatch_.initInterpolate
         (
@@ -741,6 +760,9 @@ void Foam::cyclicACMIFvPatchField<Type>::updateInterfaceMatrix
                 recvRequests_,
                 recvBufs_
             );
+
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
     }
     else
     {

src/finiteVolume/fields/fvPatchFields/constraint/cyclicAMI/cyclicAMIFvPatchField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -446,14 +446,27 @@ void Foam::cyclicAMIFvPatchField<Type>::initEvaluate
 
         // Start sending
 
-        // By-pass polyPatch to get nbrId. Instead use cyclicAMIFvPatch virtual
-        // neighbPatch()
+        // Bypass polyPatch to get nbrId.
+        // - use cyclicACMIFvPatch::neighbPatch() virtual instead
         const cyclicAMIFvPatch& neighbPatch = cyclicAMIPatch_.neighbPatch();
         const labelUList& nbrFaceCells = neighbPatch.faceCells();
         const Field<Type> pnf(this->primitiveField(), nbrFaceCells);
 
         const cyclicAMIPolyPatch& cpp = cyclicAMIPatch_.cyclicAMIPatch();
 
+        // Assert that all receives are known to have finished
+        if (!recvRequests_.empty())
+        {
+            FatalErrorInFunction
+                << "Outstanding recv request(s) on patch "
+                << cyclicAMIPatch_.name()
+                << " field " << this->internalField().name()
+                << abort(FatalError);
+        }
+
+        // Assume that sends are also OK
+        sendRequests_.clear();
+
         cpp.initInterpolate
         (
             pnf,
@@ -516,6 +529,10 @@ void Foam::cyclicAMIFvPatchField<Type>::evaluate
                 defaultValues
             ).ptr()
         );
+
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
+
         auto& patchNeighbourField = patchNeighbourFieldPtr_.ref();
 
         if (doTransform())
@@ -563,6 +580,19 @@ void Foam::cyclicAMIFvPatchField<Type>::initInterfaceMatrixUpdate
 
         const cyclicAMIPolyPatch& cpp = cyclicAMIPatch_.cyclicAMIPatch();
 
+        // Assert that all receives are known to have finished
+        if (!recvRequests_.empty())
+        {
+            FatalErrorInFunction
+                << "Outstanding recv request(s) on patch "
+                << cyclicAMIPatch_.name()
+                << " field " << this->internalField().name()
+                << abort(FatalError);
+        }
+
+        // Assume that sends are also OK
+        sendRequests_.clear();
+
         cpp.initInterpolate
         (
             pnf,
@@ -624,6 +654,9 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
                 scalarRecvBufs_,
                 defaultValues
             );
+
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
     }
     else
     {
@@ -682,6 +715,19 @@ void Foam::cyclicAMIFvPatchField<Type>::initInterfaceMatrixUpdate
 
         const cyclicAMIPolyPatch& cpp = cyclicAMIPatch_.cyclicAMIPatch();
 
+        // Assert that all receives are known to have finished
+        if (!recvRequests_.empty())
+        {
+            FatalErrorInFunction
+                << "Outstanding recv request(s) on patch "
+                << cyclicAMIPatch_.name()
+                << " field " << this->internalField().name()
+                << abort(FatalError);
+        }
+
+        // Assume that sends are also OK
+        sendRequests_.clear();
+
         cpp.initInterpolate
         (
             pnf,
@@ -742,6 +788,9 @@ void Foam::cyclicAMIFvPatchField<Type>::updateInterfaceMatrix
                 recvBufs_,
                 defaultValues
             );
+
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
     }
     else
     {

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicACMIGAMGInterfaceField/cyclicACMIGAMGInterfaceField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013 OpenFOAM Foundation
-    Copyright (C) 2019,2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -217,6 +217,18 @@ void Foam::cyclicACMIGAMGInterfaceField::initInterfaceMatrixUpdate
           : AMI.srcMap()
         );
 
+        // Assert that all receives are known to have finished
+        if (!recvRequests_.empty())
+        {
+            FatalErrorInFunction
+                << "Outstanding recv request(s) on patch "
+                << cyclicACMIInterface_.index()
+                << abort(FatalError);
+        }
+
+        // Assume that sends are also OK
+        sendRequests_.clear();
+
         // Insert send/receive requests (non-blocking). See e.g.
         // cyclicAMIPolyPatchTemplates.C
         const label oldWarnComm = UPstream::commWarn(AMI.comm());
@@ -276,6 +288,9 @@ void Foam::cyclicACMIGAMGInterfaceField::updateInterfaceMatrix
         solveScalarField work;
         map.receive(recvRequests_, scalarRecvBufs_, work);
 
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
+
         solveScalarField pnf(faceCells.size(), Zero);
         AMI.weightedSum
         (

src/meshTools/AMIInterpolation/GAMG/interfaceFields/cyclicAMIGAMGInterfaceField/cyclicAMIGAMGInterfaceField.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2019,2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -218,6 +218,18 @@ void Foam::cyclicAMIGAMGInterfaceField::initInterfaceMatrixUpdate
           : AMI.srcMap()
         );
 
+        // Assert that all receives are known to have finished
+        if (!recvRequests_.empty())
+        {
+            FatalErrorInFunction
+                << "Outstanding recv request(s) on patch "
+                << cyclicAMIInterface_.index()
+                << abort(FatalError);
+        }
+
+        // Assume that sends are also OK
+        sendRequests_.clear();
+
         // Insert send/receive requests (non-blocking). See e.g.
         // cyclicAMIPolyPatchTemplates.C
         const label oldWarnComm = UPstream::commWarn(AMI.comm());
@@ -290,6 +302,9 @@ void Foam::cyclicAMIGAMGInterfaceField::updateInterfaceMatrix
         solveScalarField work;
         map.receive(recvRequests_, scalarRecvBufs_, work);
 
+        // Receive requests all handled by last function call
+        recvRequests_.clear();
+
         solveScalarField pnf(faceCells.size(), Zero);
         AMI.weightedSum
         (

src/meshTools/triangulatedPatch/triangulatedPatch.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,7 +26,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "triangulatedPatch.H"
-#include "triPointRef.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -50,6 +49,8 @@ bool Foam::triangulatedPatch::randomPoint
 
     // Find corresponding decomposed face triangle
     // Note: triWght_ is sized nTri+1 (zero added at start)
+    //
+    // TBD: binary search with findLower(triWght_, c) ??
     label trii = 0;
     for (label i = 0; i < triWght_.size() - 1; ++i)
     {
@@ -62,11 +63,9 @@ bool Foam::triangulatedPatch::randomPoint
 
     // Find random point in triangle
     const pointField& points = patch_.points();
-    const face& tf = triFace_[trii];
-    const triPointRef tri(points[tf[0]], points[tf[1]], points[tf[2]]);
 
-    result = tri.randomPoint(rnd);
-    facei = triToFace_[trii];
+    result = triFace_[trii].tri(points).randomPoint(rnd);
+    facei = triFace_[trii].index();
     celli = patch_.faceCells()[facei];
 
     if (perturbTol_ > 0)
@@ -97,7 +96,6 @@ Foam::triangulatedPatch::triangulatedPatch
     patch_(patch),
     perturbTol_(perturbTol),
     triFace_(),
-    triToFace_(),
     triWght_()
 {
     update();
@@ -115,70 +113,95 @@ Foam::triangulatedPatch::triangulatedPatch
 {}
 
 
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+void Foam::triangulatedPatch::triangulate
+(
+    const polyPatch& pp,
+    List<labelledTri>& tris
+)
+{
+    const pointField& points = pp.points();
+
+    // Triangulate the patch faces and create addressing
+    label nTris = 0;
+    for (const face& f : pp)
+    {
+        nTris += f.nTriangles();
+    }
+
+    DynamicList<labelledTri> dynTris(nTris);
+    DynamicList<face> tfaces(8);  // work array
+
+    label facei = 0;
+    for (const face& f : pp)
+    {
+        tfaces.clear();
+        f.triangles(points, tfaces);
+
+        for (const auto& t : tfaces)
+        {
+            dynTris.emplace_back(t[0], t[1], t[2], facei);
+        }
+        ++facei;
+    }
+
+    tris.transfer(dynTris);
+}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::triangulatedPatch::update()
 {
+    triFace_.clear();
+    triWght_.clear();
+
+    triangulate(patch_, triFace_);
+
     const pointField& points = patch_.points();
 
-    // Triangulate the patch faces and create addressing
-    DynamicList<label> triToFace(2*patch_.size());
-    DynamicList<face> triFace(2*patch_.size());
-    DynamicList<scalar> triWght(2*patch_.size());
-    DynamicList<face> tris(8);
+    const label myProci = UPstream::myProcNo();
+    const label numProc = UPstream::nProcs();
+
+    // Calculate the cumulative triangle weights
+    triWght_.resize_nocopy(triFace_.size()+1);
+
+    auto iter = triWght_.begin();
 
     // Set zero value at the start of the tri area/weight list
-    triWght.push_back(0);
+    scalar patchArea = 0;
+    *iter = patchArea;
+    ++iter;
 
-    forAll(patch_, facei)
+    // Calculate cumulative and total area (processor-local at this point)
+    for (const auto& t : triFace_)
     {
-        const face& f = patch_[facei];
+        patchArea += t.mag(points);
 
-        tris.clear();
-        f.triangles(points, tris);
-
-        for (const auto& t : tris)
-        {
-            triToFace.push_back(facei);
-            triFace.push_back(t);
-            triWght.push_back(t.mag(points));
-        }
+        *iter = patchArea;
+        ++iter;
     }
 
-    scalarList procSumWght(Pstream::nProcs()+1, Zero);
-    procSumWght[Pstream::myProcNo()+1] = sum(triWght);
+    // FIXME: use allGatherList of subslice
+    scalarList procSumWght(numProc+1, Foam::zero{});
+    procSumWght[myProci+1] = patchArea;
     Pstream::listCombineReduce(procSumWght, maxEqOp<scalar>());
 
+    // Convert to cumulative
     for (label i = 1; i < procSumWght.size(); ++i)
     {
-        // Convert to cumulative
         procSumWght[i] += procSumWght[i-1];
     }
 
-    const scalar offset = procSumWght[Pstream::myProcNo()];
-    forAll(triWght, i)
+    const scalar offset = procSumWght[myProci];
+    const scalar totalArea = procSumWght.back();
+
+    // Apply processor offset and normalise - for a global 0-1 interval
+    for (scalar& w : triWght_)
     {
-        if (i)
-        {
-            // Convert to cumulative
-            triWght[i] += triWght[i-1];
-        }
-
-        // Apply processor offset
-        triWght[i] += offset;
+        w = (w + offset) / totalArea;
     }
-
-    // Normalise
-    const scalar sumWght = procSumWght.back();
-    for (scalar& w : triWght)
-    {
-        w /= sumWght;
-    }
-
-    // Transfer to persistent storage
-    triFace_.transfer(triFace);
-    triToFace_.transfer(triToFace);
-    triWght_.transfer(triWght);
 }
 
 
@@ -190,6 +213,14 @@ bool Foam::triangulatedPatch::randomLocalPoint
     label& celli
 ) const
 {
+    if (triWght_.empty())
+    {
+        result = point::min;
+        facei = -1;
+        celli = -1;
+        return false;
+    }
+
     const scalar c = rnd.position<scalar>(triWght_.front(), triWght_.back());
 
     return randomPoint(rnd, c, result, facei, celli);
@@ -204,20 +235,23 @@ bool Foam::triangulatedPatch::randomGlobalPoint
     label& celli
 ) const
 {
-    boolList valid(UPstream::nProcs(), false);
-    valid[UPstream::myProcNo()] = randomLocalPoint(rnd, result, facei, celli);
-    UPstream::listGatherValues(valid);
-
-    forAll(valid, proci)
+    if (UPstream::parRun())
     {
-        // Choose first valid processor
-        if (valid[proci])
-        {
-            return (proci == UPstream::myProcNo());
-        }
-    }
+        const scalar c = rnd.sample01<scalar>();
+        const bool ok = randomPoint(rnd, c, result, facei, celli);
 
-    return false;
+        boolList valid(UPstream::listGatherValues(ok));
+
+        // Select the first valid processor
+        label proci = valid.find(true);
+        Pstream::broadcast(proci);
+
+        return (proci == UPstream::myProcNo());
+    }
+    else
+    {
+        return randomLocalPoint(rnd, result, facei, celli);
+    }
 }
 
 

src/meshTools/triangulatedPatch/triangulatedPatch.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2023 OpenCFD Ltd.
+    Copyright (C) 2023-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,6 +35,7 @@ Description
 #define Foam_triangulatedPatch_H
 
 #include "polyMesh.H"
+#include "labelledTri.H"
 #include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -55,13 +56,12 @@ class triangulatedPatch
         //- Perturbation tolerance to move the point towards the cell centre
         bool perturbTol_;
 
-        //- Face triangles
-        faceList triFace_;
+        //- The polyPatch faces as triangles, the index of each corresponds
+        //- to the undecomposed patch face index.
+        List<labelledTri> triFace_;
 
-        //- Triangle to patch face addressing
-        labelList triToFace_;
-
-        //- Triangle weights
+        //- The cumulative triangle area per triangle face,
+        //- globally normalised as a 0-1 interval.
         scalarList triWght_;
 
 
@@ -81,6 +81,11 @@ class triangulatedPatch
         ) const;
 
 
+    // Static Member Functions
+
+        //- Triangulate the patch faces and create addressing
+        static void triangulate(const polyPatch& pp, List<labelledTri>& tris);
+
 public:
 
     //- Constructors

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_;