ENH: pointMesh: support for subsetMesh

- surfaceWriter TryNew() factory methods for more failure tolerant
  handling

- reduce communication for sampledSurfaces.
  Track non-empty surfaces as bool, only updated on change
  (expire/update).

CONTRIBUTORS.md

@@ -9,6 +9,7 @@ It is likely incomplete...
 - Yu Ankun
 - Tetsuo Aoyagi
 - Akira Azami
+- Pete Bachant
 - William Bainbridge
 - Gabriel Barajas
 - Kutalmis Bercin

applications/test/PtrList/Test-PtrList.C

@@ -412,7 +412,7 @@ int main(int argc, char *argv[])
     }
 
     {
-        Info<< "range-for of list (" << list1.count() << '/'
+        Info<< "range-for of list (" << list1.count_nonnull() << '/'
             << list1.size() << ") non-null entries" << nl
             << "(" << nl;
         for (const auto& item : list1)

applications/test/syncTools/missing_cell/Allclean

@@ -0,0 +1,11 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
+#------------------------------------------------------------------------------
+
+cleanCase
+
+# Remove surface and features
+rm -rf constant/triSurface
+
+#------------------------------------------------------------------------------

applications/test/syncTools/missing_cell/Allrun

@@ -0,0 +1,27 @@
+#!/bin/sh
+cd "${0%/*}" || exit                                # Run from this directory
+. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
+#------------------------------------------------------------------------------
+
+#- Generate 2x2x1 cells
+runApplication blockMesh
+
+#- Remove cell0
+runApplication topoSet
+runApplication subsetMesh c0 -patch exposed0 -overwrite
+
+#- Put exposed faces (2) into separate patches
+runApplication -s face topoSet
+runApplication createPatch -overwrite
+
+#- Decompose - creates one processor without any faces in patches
+runApplication decomposePar
+
+#- Extract inter-patch points. Should include processor that does not
+#- have faces on patch ...
+mkdir -p constant/triSurface
+runParallel surfaceMeshExtract \
+    -patches '(exposed0 exposed1)' -featureAngle 180 \
+    constant/triSurface/blockMesh.obj
+
+#------------------------------------------------------------------------------

applications/test/syncTools/missing_cell/README.txt

@@ -0,0 +1,7 @@
+- 2x2x1 mesh
+- remove one cell, exposing two faces
+- move exposed faces into two patches
+- decompose onto 3
+- run surfaceMeshExtract -featureAngle 180
+- should also mark points on the processor that has no
+  faces but is coupled

applications/test/syncTools/missing_cell/constant/transportProperties

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2312                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+nu              0.01;
+
+
+// ************************************************************************* //

applications/test/syncTools/missing_cell/system/blockMeshDict

@@ -0,0 +1,88 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2312                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scale   1;
+
+vertices
+(
+    //- Single block
+    (0  0   0)
+    (2  0   0)
+    (2  2   0)
+    (0  2   0)
+    (0  0   2)
+    (2  0   2)
+    (2  2   2)
+    (0  2   2)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (2 2 1) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    topWall
+    {
+        type wall;
+        faces
+        (
+            (3 7 6 2)
+        );
+    }
+    bottomWall
+    {
+        type wall;
+        faces
+        (
+            (1 5 4 0)
+        );
+    }
+    fixedWalls
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+            (2 6 5 1)
+        );
+    }
+    frontAndBack
+    {
+        type patch;
+        faces
+        (
+            (0 3 2 1)
+            (4 5 6 7)
+        );
+    }
+    exposed0
+    {
+        type patch;
+        faces   ();
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

applications/test/syncTools/missing_cell/system/controlDict

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2312                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     icoFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         0.5;
+
+deltaT          0.005;
+
+writeControl    timeStep;
+
+writeInterval   20;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  16;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+// ************************************************************************* //

applications/test/syncTools/missing_cell/system/createPatchDict

@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2312                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      createPatchDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+pointSync false;
+
+// Patches to create.
+patches
+(
+    // Example of creating mapped patches using geometric matching
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    {
+        // Name of new patch
+        name exposed1;
+
+        // Dictionary to construct new patch from
+        patchInfo
+        {
+            type        patch;
+        }
+
+        // How to select the faces:
+        //  - set : specify faceSet in 'set'
+        //  - patches : specify names in 'patches'
+        //  - autoPatch : attempts automatic patching of the specified
+        //                candidates in 'patches'.
+        constructFrom   set;
+        set             exposed0;
+    }
+);
+
+// ************************************************************************* //

applications/test/syncTools/missing_cell/system/decomposeParDict

@@ -0,0 +1,24 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2312                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    note        "mesh decomposition control dictionary";
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- The total number of domains (mandatory)
+numberOfSubdomains  3;
+
+//- The decomposition method (mandatory)
+method              scotch;
+
+// ************************************************************************* //

applications/test/syncTools/missing_cell/system/fvSchemes

@@ -0,0 +1,51 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2312                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+    default         Euler;
+}
+
+gradSchemes
+{
+    default         Gauss linear;
+    grad(p)         Gauss linear;
+}
+
+divSchemes
+{
+    default         none;
+    div(phi,U)      Gauss linear;
+}
+
+laplacianSchemes
+{
+    default         Gauss linear orthogonal;
+}
+
+interpolationSchemes
+{
+    default         linear;
+}
+
+snGradSchemes
+{
+    default         orthogonal;
+}
+
+
+// ************************************************************************* //

applications/test/syncTools/missing_cell/system/fvSolution

@@ -0,0 +1,52 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2312                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+    p
+    {
+        solver          PCG;
+        preconditioner  DIC;
+        tolerance       1e-06;
+        relTol          0.05;
+    }
+
+    pFinal
+    {
+        $p;
+        relTol          0;
+    }
+
+    U
+    {
+        solver          smoothSolver;
+        smoother        symGaussSeidel;
+        tolerance       1e-05;
+        relTol          0;
+    }
+}
+
+PISO
+{
+    nCorrectors     2;
+    nNonOrthogonalCorrectors 0;
+    pRefCell        0;
+    pRefValue       0;
+}
+
+
+// ************************************************************************* //

applications/test/syncTools/missing_cell/system/topoSetDict

@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2312                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      topoSetDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+actions
+(
+    {
+        name    c0;
+        type    cellSet;
+        action  new;
+        source  labelToCell;
+        value   (0);
+    }
+    {
+        name    c0;
+        type    cellSet;
+        action  invert;
+    }
+
+    {
+        name    exposed0;
+        type    faceSet;
+        action  new;
+        source  patchToFace;
+        patch   exposed0;
+    }
+    {
+        name    exposed0;
+        type    faceSet;
+        action  subset;
+        source  boxToFace;
+        box     (-100 1 -100)(100 100 100);
+    }
+);
+
+// ************************************************************************* //

applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.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.
@@ -220,7 +220,7 @@ PtrList<FieldType> subsetFields
     const pointMesh& pMesh
 )
 {
-    const fvMesh& baseMesh = subsetter.baseMesh();
+    //const fvMesh& baseMesh = subsetter.baseMesh();
 
     const UPtrList<const IOobject> fieldObjects
     (
@@ -247,8 +247,8 @@ PtrList<FieldType> subsetFields
             IOobject
             (
                 io.name(),
-                baseMesh.time().timeName(),
-                baseMesh,
+                pMesh.thisDb().time().timeName(),
+                pMesh.thisDb(),
                 IOobjectOption::MUST_READ,
                 IOobjectOption::NO_WRITE,
                 IOobjectOption::NO_REGISTER
@@ -382,6 +382,8 @@ int main(int argc, char *argv[])
     #include "createTime.H"
 
     #include "createNamedMesh.H"
+    // Make sure pointMesh gets constructed/read as well
+    (void)pointMesh::New(mesh, IOobject::READ_IF_PRESENT);
 
     // arg[1] = word (cellSet) or wordRes (cellZone)
     // const word selectionName = args[1];
@@ -583,7 +585,7 @@ int main(int argc, char *argv[])
     // Read point fields and subset
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    const pointMesh& pMesh = pointMesh::New(mesh);
+    const pointMesh& pMesh = pointMesh::New(mesh, IOobject::READ_IF_PRESENT);
 
     #undef  createSubsetFields
     #define createSubsetFields(FieldType, Variable)             \
@@ -663,6 +665,18 @@ int main(int argc, char *argv[])
     subsetter.subMesh().write();
     processorMeshes::removeFiles(subsetter.subMesh());
 
+    auto* subPointMeshPtr =
+        subsetter.subMesh().thisDb().findObject<pointMesh>
+        (
+            pointMesh::typeName
+        );
+    if (subPointMeshPtr)
+    {
+        pointMesh& subPointMesh = const_cast<pointMesh&>(*subPointMeshPtr);
+        subPointMesh.setInstance(subsetter.subMesh().facesInstance());
+        subPointMesh.write();
+    }
+
 
     // Volume fields
     for (const auto& fld : vScalarFlds)     { fld.write(); }

applications/utilities/parallelProcessing/decomposePar/decomposePar.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -661,6 +661,9 @@ int main(int argc, char *argv[])
             ),
             decompDictFile
         );
+        // Make sure pointMesh gets read as well
+        (void)pointMesh::New(mesh, IOobject::READ_IF_PRESENT);
+
 
         // Decompose the mesh
         if (!decomposeFieldsOnly)
@@ -780,6 +783,7 @@ int main(int argc, char *argv[])
             PtrList<labelIOList> cellProcAddressingList(mesh.nProcs());
             PtrList<labelIOList> boundaryProcAddressingList(mesh.nProcs());
             PtrList<labelIOList> pointProcAddressingList(mesh.nProcs());
+            PtrList<labelIOList> pointBoundaryProcAddressingList(mesh.nProcs());
 
             PtrList<fvFieldDecomposer> fieldDecomposerList(mesh.nProcs());
             PtrList<pointFieldDecomposer> pointFieldDecomposerList
@@ -850,7 +854,10 @@ int main(int argc, char *argv[])
 
                 // Point fields
                 // ~~~~~~~~~~~~
-                const pointMesh& pMesh = pointMesh::New(mesh);
+
+                // Read decomposed pointMesh
+                const pointMesh& pMesh =
+                    pointMesh::New(mesh, IOobject::READ_IF_PRESENT);
 
                 pointFieldDecomposer::fieldsCache pointFieldCache;
 
@@ -1119,7 +1126,34 @@ int main(int argc, char *argv[])
                             pointProcAddressingList
                         );
 
-                        const pointMesh& procPMesh = pointMesh::New(procMesh);
+                        const pointMesh& procPMesh =
+                            pointMesh::New(procMesh, IOobject::READ_IF_PRESENT);
+
+                        if (!pointBoundaryProcAddressingList.set(proci))
+                        {
+                            pointBoundaryProcAddressingList.set
+                            (
+                                proci,
+                                autoPtr<labelIOList>::New
+                                (
+                                    IOobject
+                                    (
+                                        "boundaryProcAddressing",
+                                        procMesh.facesInstance(),
+                                        polyMesh::meshSubDir
+                                       /pointMesh::meshSubDir,
+                                        procPMesh.thisDb(),
+                                        IOobject::READ_IF_PRESENT,
+                                        IOobject::NO_WRITE,
+                                        IOobject::NO_REGISTER
+                                    ),
+                                    boundaryProcAddressing
+                                )
+                            );
+                        }
+                        const auto& pointBoundaryProcAddressing =
+                            pointBoundaryProcAddressingList[proci];
+
 
                         if (!pointFieldDecomposerList.set(proci))
                         {
@@ -1131,7 +1165,7 @@ int main(int argc, char *argv[])
                                     pMesh,
                                     procPMesh,
                                     pointProcAddressing,
-                                    boundaryProcAddressing
+                                    pointBoundaryProcAddressing
                                 )
                             );
                         }
@@ -1143,6 +1177,12 @@ int main(int argc, char *argv[])
 
                         if (times.size() == 1)
                         {
+                            // Early deletion
+                            pointBoundaryProcAddressingList.set
+                            (
+                                proci,
+                                nullptr
+                            );
                             pointProcAddressingList.set(proci, nullptr);
                             pointFieldDecomposerList.set(proci, nullptr);
                         }

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -44,6 +44,12 @@ License
 #include "decompositionModel.H"
 #include "hexRef8Data.H"
 
+// For handling pointMeshes with additional patches
+#include "pointMesh.H"
+#include "meshPointPatch.H"
+#include "processorPointPatch.H"
+#include "DynamicField.H"
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::domainDecomposition::mark
@@ -740,6 +746,101 @@ bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
 
         procMesh.write();
 
+        // Add pointMesh if it was available
+        const auto* pMeshPtr =
+            thisDb().cfindObject<pointMesh>(pointMesh::typeName);
+        if (pMeshPtr)
+        {
+            const auto& pMesh = *pMeshPtr;
+            const auto& pMeshBoundary = pMesh.boundary();
+
+
+            // 1. Generate pointBoundaryMesh from polyBoundaryMesh (so ignoring
+            //    any additional patches
+            const auto& procPointMesh = pointMesh::New(procMesh);
+
+            pointBoundaryMesh& procBoundary =
+                const_cast<pointBoundaryMesh&>(procPointMesh.boundary());
+
+
+            // 2. Explicitly add subsetted meshPointPatches
+            forAll(pMeshBoundary, patchi)
+            {
+                const auto* mppPtr = isA<meshPointPatch>(pMeshBoundary[patchi]);
+                if (mppPtr && (procBoundary.findPatchID(mppPtr->name()) == -1))
+                {
+                    const auto& mpp = *mppPtr;
+
+                    DynamicList<label> procMeshPoints(mpp.size());
+                    DynamicField<vector> procNormals(mpp.size());
+                    forAll(mpp.meshPoints(), i)
+                    {
+                        const label pointi = mpp.meshPoints()[i];
+                        const label procPointi = pointLookup[pointi];
+                        if (procPointi != -1)
+                        {
+                            procMeshPoints.append(procPointi);
+                            procNormals.append(mpp.pointNormals()[i]);
+                        }
+                    }
+
+                    procBoundary.push_back
+                    (
+                        new meshPointPatch
+                        (
+                            mpp.name(),
+                            procMeshPoints,
+                            procNormals,
+                            procBoundary.size(),
+                            procBoundary,
+                            meshPointPatch::typeName
+                        )
+                    );
+                }
+            }
+
+            // 3. Shuffle new patches before any processor patches
+            labelList oldToNew(procBoundary.size());
+            label newPatchi = 0;
+            forAll(procBoundary, patchi)
+            {
+                if (!isA<processorPointPatch>(procBoundary[patchi]))
+                {
+                    oldToNew[patchi] = newPatchi;
+                    newPatchi++;
+                }
+            }
+
+            // decomposed-to-undecomposed patch numbering
+            labelList boundaryProcAddressing(identity(newPatchi));
+            boundaryProcAddressing.setSize(procBoundary.size(), -1);
+
+            forAll(procBoundary, patchi)
+            {
+                if (isA<processorPointPatch>(procBoundary[patchi]))
+                {
+                    oldToNew[patchi] = newPatchi++;
+                }
+            }
+            procBoundary.reorder(oldToNew, true);
+
+            // Write pointMesh/boundary
+            procBoundary.write();
+
+            // Write pointMesh/boundaryProcAddressing
+            IOobject ioAddr
+            (
+                "boundaryProcAddressing",
+                procMesh.facesInstance(),
+                polyMesh::meshSubDir/pointMesh::meshSubDir,
+                procPointMesh.thisDb(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                IOobject::NO_REGISTER
+            );
+            IOListRef<label>(ioAddr, boundaryProcAddressing).write();
+        }
+
         // Write points if pointsInstance differing from facesInstance
         if (facesInstancePointsPtr_)
         {

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -411,24 +411,18 @@ int main(int argc, char *argv[])
             {
                 Info<< "Reconstructing point fields" << nl << endl;
 
-                const pointMesh& pMesh = pointMesh::New(mesh);
-                PtrList<pointMesh> pMeshes(procMeshes.meshes().size());
-
-                forAll(pMeshes, proci)
-                {
-                    pMeshes.set
-                    (
-                        proci,
-                        new pointMesh(procMeshes.meshes()[proci])
-                    );
-                }
+                const pointMesh& pMesh = pointMesh::New
+                (
+                    mesh,
+                    IOobject::READ_IF_PRESENT
+                );
 
                 pointFieldReconstructor reconstructor
                 (
                     pMesh,
-                    pMeshes,
+                    procMeshes.pointMeshes(),
                     procMeshes.pointProcAddressing(),
-                    procMeshes.boundaryProcAddressing()
+                    procMeshes.pointMeshBoundaryProcAddressing()
                 );
 
                 reconstructor.reconstructAllFields(objects, selectedFields);

applications/utilities/surface/surfaceMeshExtract/surfaceMeshExtract.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2023 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,10 +34,13 @@ Description
     Extract patch or faceZone surfaces from a polyMesh.
     Depending on output surface format triangulates faces.
 
-    Region numbers on faces no guaranteed to be the same as the patch indices.
+    Region numbers on faces not guaranteed to be the same as the patch indices.
 
     Optionally only extracts named patches.
 
+    Optionally filters out points on feature-edges and generates pointPatches
+    for these - written to pointMesh/boundary.
+
     If run in parallel, processor patches get filtered out by default and
     the mesh is merged (based on topology).
 
@@ -48,6 +51,7 @@ Description
 #include "argList.H"
 #include "Time.H"
 #include "polyMesh.H"
+#include "pointMesh.H"
 #include "emptyPolyPatch.H"
 #include "processorPolyPatch.H"
 #include "ListListOps.H"
@@ -55,6 +59,11 @@ Description
 #include "globalMeshData.H"
 #include "globalIndex.H"
 #include "timeSelector.H"
+#include "meshPointPatch.H"
+#include "unitConversion.H"
+#include "dummyTransform.H"
+#include "syncTools.H"
+#include "processorPointPatch.H"
 
 using namespace Foam;
 
@@ -106,6 +115,551 @@ labelList getSelectedPatches
 }
 
 
+label addMeshPointPatches
+(
+    const polyMesh& mesh,
+    const pointMesh& pMesh,
+
+    const uindirectPrimitivePatch& allBoundary,
+    const labelUList& faceToZone,
+    const surfZoneIdentifierList& surfZones,
+
+    const scalar edgeFeatureAngle,
+    const scalar pointFeatureAngle,
+    const bool verbose = true
+)
+{
+    const pointBoundaryMesh& pointBm = pMesh.boundary();
+    const label nPointPatches = pointBm.size();
+
+    const globalMeshData& globalData = mesh.globalData();
+    const indirectPrimitivePatch& cpp = globalData.coupledPatch();
+    const mapDistribute& map = globalData.globalEdgeSlavesMap();
+    const auto& mp = allBoundary.meshPoints();
+
+
+    const vector nullVector(vector::uniform(0));
+    const auto assignNonNull = [&](vector& x, const vector& y)
+    {
+        if (x == nullVector && y != nullVector)
+        {
+            x = y;
+        }
+    };
+
+    // Calculate parallel-consistent point normals (as unweighted average
+    // of faceNormals). Note: only valid on patch points, not on mesh points
+    // that are coupled to these.
+    const pointField pointNormals
+    (
+        PatchTools::pointNormals
+        (
+            mesh,
+            allBoundary
+        )
+    );
+    // Expand to all coupled points
+    pointField meshPointNormals(mesh.nPoints(), nullVector);
+    UIndirectList<vector>(meshPointNormals, mp) = pointNormals;
+    syncTools::syncPointList
+    (
+        mesh,
+        meshPointNormals,
+        assignNonNull,
+        nullVector
+    );
+
+
+    // Find correspondence between allBoundary and coupled edges
+    labelList allEdges;
+    labelList coupledEdges;
+    bitSet sameEdgeOrientation;
+    PatchTools::matchEdges
+    (
+        allBoundary,
+        cpp,
+        allEdges,
+        coupledEdges,
+        sameEdgeOrientation
+    );
+
+
+
+    // To construct the patches we need to know per edge
+    //  - patch on either side (if topological feature edge)
+    //  - faceNormal on either side (if feature angle)
+    // We need to know per point
+    //  - patches on all connected faces
+    //  - faceNormals on all connected faces? And compare to average?
+    //    or edge normals on all connected edges
+
+
+    typedef Tuple2<label, vector> PN;
+    const PN nullPN(-1, nullVector);
+
+
+
+    // Point-based analysis
+    // ~~~~~~~~~~~~~~~~~~~~
+
+
+    // Collect per (mesh)point the zones (1, 2 or >2). Note: per mesh to
+    // make it easier to sync. See edge-based code below where we explicitly
+    // have to transfer from patch-edge to mesh-point etc. Note sure which one
+    // fits better....
+    labelPairList pointToZones(mesh.nPoints(), labelPair(-1, -1));
+    {
+        // Combine zones.
+        const auto combineZones = [&](labelPair& x, const labelPair& y)
+        {
+            if (x == labelPair(-2, -2))
+            {
+                // Already marked
+            }
+            else if (y == labelPair(-2, -2))
+            {
+                x = y;
+            }
+            else
+            {
+                // Find first free slot
+                if (x[0] == -1)
+                {
+                    if (y[0] != -1)
+                    {
+                        x[0] = y[0];
+                    }
+                    else
+                    {
+                        x[0] = y[1];
+                    }
+                }
+                else if (x[1] == -1)
+                {
+                    if (y[0] != -1 && y[0] != x[0])
+                    {
+                        x[1] = y[0];
+                    }
+                    else if (y[1] != -1 && y[1] != x[1])
+                    {
+                        x[1] = y[1];
+                    }
+                }
+                else
+                {
+                    // Both x slots filled. See if y adds a 3rd element
+                    if (y[0] != -1 && y[0] != x[0] && y[0] != x[1])
+                    {
+                        x = labelPair(-2, -2);
+                    }
+                    else if (y[1] != -1 && y[1] != x[0] && y[1] != x[1])
+                    {
+                        x = labelPair(-2, -2);
+                    }
+                }
+            }
+        };
+
+
+        forAll(allBoundary, facei)
+        {
+            const auto& f = allBoundary[facei];
+            const label zonei = faceToZone[facei];
+            for (const label pointi : f)
+            {
+                auto& pZones = pointToZones[pointi];
+
+                if (pZones != labelPair(-2, -2) && !pZones.contains(zonei))
+                {
+                    if (pZones.first() == -1)
+                    {
+                        pZones.first() = zonei;
+                    }
+                    else if (pZones.second() == -1)
+                    {
+                        pZones.second() = zonei;
+                    }
+                    else
+                    {
+                        // Mark as >2 zones
+                        pZones = labelPair(-2, -2);
+                    }
+                }
+            }
+        }
+
+        syncTools::syncPointList
+        (
+            mesh,
+            pointToZones,
+            combineZones,
+            labelPair(-1, -1),
+            dummyTransform()
+        );
+    }
+
+
+
+
+    // Edge-based analysis
+    // ~~~~~~~~~~~~~~~~~~~~
+
+    // 1. Local analysis
+
+    List<Pair<PN>> allEdgeToFaces
+    (
+        allBoundary.nEdges(),
+        Pair<PN>(nullPN, nullPN)
+    );
+    {
+        const auto& edgeFaces = allBoundary.edgeFaces();
+        const auto& faceNormals = allBoundary.faceNormals();
+
+        forAll(edgeFaces, edgei)
+        {
+            const auto& eFaces = edgeFaces[edgei];
+            const vector& n0 = faceNormals[eFaces[0]];
+            const label zone0 = faceToZone[eFaces[0]];
+            if (eFaces.size() == 1)
+            {
+                allEdgeToFaces[edgei] = Pair<PN>(PN(zone0, n0), nullPN);
+            }
+            else
+            {
+                const vector& n1 = faceNormals[eFaces[1]];
+                const label zone1 = faceToZone[eFaces[1]];
+                allEdgeToFaces[edgei] = Pair<PN>
+                (
+                    PN(zone0, n0),
+                    PN(zone1, n1)
+                );
+            }
+        }
+    }
+
+
+    // 2. Sync across coupled patches
+
+    {
+        // Combine pair of normals
+        const auto vectorPairMax = [&](Pair<PN>& x, const Pair<PN>& y)
+        {
+            if (x[0] == nullPN)
+            {
+                if (y[0] != nullPN)
+                {
+                    x[0] = y[0];
+                }
+                else
+                {
+                    x[0] = y[1];
+                }
+            }
+            else if (x[1] == nullPN)
+            {
+                if (y[0] != nullPN && y[0] != x[0])
+                {
+                    x[1] = y[0];
+                }
+                else
+                {
+                    x[1] = y[1];
+                }
+            }
+        };
+
+        List<Pair<PN>> cppEdgeData
+        (
+            map.constructSize(),
+            Pair<PN>(nullPN, nullPN)
+        );
+        UIndirectList<Pair<PN>>(cppEdgeData, coupledEdges) =
+            UIndirectList<Pair<PN>>(allEdgeToFaces, allEdges);
+
+        globalData.syncData
+        (
+            cppEdgeData,
+            globalData.globalEdgeSlaves(),
+            globalData.globalEdgeTransformedSlaves(),
+            map,
+            globalData.globalTransforms(),
+            vectorPairMax,
+            dummyTransform()
+        );
+
+        UIndirectList<Pair<PN>>(allEdgeToFaces, allEdges) =
+            UIndirectList<Pair<PN>>(cppEdgeData, coupledEdges);
+    }
+
+
+    // Now we have all the per-patch edge information
+    // - do inter-patch edges
+    // - do feature-angle edges
+    // Store on mesh points
+
+    const auto assignNonNullPN = [&](PN& x, const PN& y)
+    {
+        if (x.second() == nullVector && y.second() != nullVector)
+        {
+            x = y;
+        }
+    };
+
+
+    // Storing the normal for points that are on inter-patch edges
+    vectorField patchEdgeNormal(mesh.nPoints(), nullVector);
+    // Storing the normal for points that are on patch-internal feat edges
+    List<PN> featEdgeNormal(mesh.nPoints(), nullPN);
+    forAll(allEdgeToFaces, edgei)
+    {
+        const edge& e = allBoundary.edges()[edgei];
+        const label mp0 = mp[e[0]];
+        const label mp1 = mp[e[1]];
+
+        const Pair<PN>& facesInfo = allEdgeToFaces[edgei];
+
+        if (facesInfo[1] == nullPN)
+        {
+            // Boundary edge
+            patchEdgeNormal[mp0] = pointNormals[e[0]];
+            patchEdgeNormal[mp1] = pointNormals[e[1]];
+        }
+        else
+        {
+            if (facesInfo[0].first() != facesInfo[1].first())
+            {
+                // Inter-patch
+                patchEdgeNormal[mp0] = pointNormals[e[0]];
+                patchEdgeNormal[mp1] = pointNormals[e[1]];
+            }
+            else
+            {
+                // Same patch - check for feature angle
+
+                const vector& n0 = facesInfo[0].second();
+                const vector& n1 = facesInfo[1].second();
+
+                if ((n0 & n1) < Foam::cos(degToRad(edgeFeatureAngle)))
+                {
+                    if (patchEdgeNormal[mp0] == nullVector)
+                    {
+                        featEdgeNormal[mp0] = PN
+                        (
+                            facesInfo[0].first(),   // zone
+                            pointNormals[e[0]]
+                        );
+                    }
+                    if (patchEdgeNormal[mp1] == nullVector)
+                    {
+                        featEdgeNormal[mp1] = PN
+                        (
+                            facesInfo[0].first(),   // zone
+                            pointNormals[e[1]]
+                        );
+                    }
+                }
+            }
+        }
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        patchEdgeNormal,
+        assignNonNull,
+        nullVector
+    );
+    syncTools::syncPointList
+    (
+        mesh,
+        featEdgeNormal,
+        assignNonNullPN,
+        nullPN,
+        dummyTransform()
+    );
+
+    // Make sure that inter-patch points are not also in feature-edge
+    // points. Note: not absolutely nessecary since all inter-patch points
+    // will also be in the 'normal' facePointPatches.
+
+    DynamicList<label> multiZoneMeshPoints(allBoundary.nPoints());
+    forAll(pointToZones, pointi)
+    {
+        if (pointToZones[pointi] == labelPair(-2, -2))
+        {
+            multiZoneMeshPoints.append(pointi);
+            // Unmark as feature angle point
+            patchEdgeNormal[pointi] = nullVector;
+            featEdgeNormal[pointi] = nullPN;
+        }
+    }
+
+
+    DynamicList<label> twoZoneMeshPoints(allBoundary.nPoints());
+    forAll(patchEdgeNormal, pointi)
+    {
+        if (patchEdgeNormal[pointi] != nullVector)
+        {
+            twoZoneMeshPoints.append(pointi);
+            // Unmark as feature angle point
+            featEdgeNormal[pointi] = nullPN;
+        }
+    }
+
+
+    // Sort featEdgeNormal according to zone
+    List<List<label>> zoneToMeshPoints(surfZones.size());
+    List<vectorField> zoneToNormal(surfZones.size());
+    {
+        labelList sizes(surfZones.size(), 0);
+        forAll(featEdgeNormal, pointi)
+        {
+            const auto& pInfo = featEdgeNormal[pointi];
+            if (pInfo != nullPN)
+            {
+                const label zonei = pInfo.first();
+                sizes[zonei]++;
+            }
+        }
+        forAll(zoneToMeshPoints, zonei)
+        {
+            zoneToMeshPoints[zonei].setSize(sizes[zonei]);
+            zoneToNormal[zonei].setSize(sizes[zonei]);
+        }
+        sizes = 0;
+        forAll(featEdgeNormal, pointi)
+        {
+            const auto& pInfo = featEdgeNormal[pointi];
+            if (pInfo != nullPN)
+            {
+                const label zonei = pInfo.first();
+                const label index = sizes[zonei]++;
+                zoneToMeshPoints[zonei][index] = pointi;
+                zoneToNormal[zonei][index] = pInfo.second();
+            }
+        }
+    }
+
+    // Add patches
+    forAll(zoneToMeshPoints, zonei)
+    {
+        const label nPoints =
+            returnReduce(zoneToMeshPoints[zonei].size(), sumOp<label>());
+
+        const word patchName(surfZones[zonei].name() + "Edges");
+
+        if (nPoints && (pointBm.findPatchID(patchName) == -1))
+        {
+            const_cast<pointBoundaryMesh&>(pointBm).push_back
+            (
+                new meshPointPatch
+                (
+                    patchName,
+                    zoneToMeshPoints[zonei],
+                    zoneToNormal[zonei],
+                    pointBm.size(),
+                    pointBm,
+                    meshPointPatch::typeName
+                )
+            );
+
+            if (verbose)
+            {
+                const auto& ppp = pointBm.last();
+                Info<< "Added feature-edges pointPatch " << ppp.name()
+                    << " with " << nPoints << " points" << endl;
+            }
+        }
+    }
+
+
+    // Add inter-patch points
+
+    const word allEdgePatchName("boundaryEdges");
+    const label nPatchEdgePoints =
+        returnReduce(twoZoneMeshPoints.size(), sumOp<label>());
+    if (nPatchEdgePoints && (pointBm.findPatchID(allEdgePatchName) == -1))
+    {
+        const_cast<pointBoundaryMesh&>(pointBm).push_back
+        (
+            new meshPointPatch
+            (
+                allEdgePatchName,
+                twoZoneMeshPoints,
+                vectorField
+                (
+                    patchEdgeNormal,    // is pointNormal expanded to all mesh
+                    twoZoneMeshPoints
+                ),
+                pointBm.size(),
+                pointBm,
+                meshPointPatch::typeName
+            )
+        );
+
+        if (verbose)
+        {
+            const auto& ppp = pointBm.last();
+            Info<< "Added inter-patch pointPatch " << ppp.name()
+                << " with " << nPatchEdgePoints << " points" << endl;
+        }
+    }
+
+
+    const word allPointPatchName("boundaryPoints");
+    const label nMultiPoints =
+        returnReduce(multiZoneMeshPoints.size(), sumOp<label>());
+    if (nMultiPoints && (pointBm.findPatchID(allPointPatchName) == -1))
+    {
+        const_cast<pointBoundaryMesh&>(pointBm).push_back
+        (
+            new meshPointPatch
+            (
+                allPointPatchName,
+                multiZoneMeshPoints,
+                vectorField
+                (
+                    meshPointNormals,   // is pointNormal expanded to all mesh
+                    multiZoneMeshPoints
+                ),
+                pointBm.size(),
+                pointBm,
+                meshPointPatch::typeName
+            )
+        );
+
+        if (verbose)
+        {
+            const auto& ppp = pointBm.last();
+            Info<< "Added multi-patch pointPatch " << ppp.name()
+                << " with " << nMultiPoints << " points" << endl;
+        }
+    }
+
+
+    // Shuffle into order
+    labelList oldToNew(pointBm.size());
+    label newPatchi = 0;
+    forAll(pointBm, patchi)
+    {
+        if (!isA<processorPointPatch>(pointBm[patchi]))
+        {
+            oldToNew[patchi] = newPatchi++;
+        }
+    }
+    forAll(pointBm, patchi)
+    {
+        if (isA<processorPointPatch>(pointBm[patchi]))
+        {
+            oldToNew[patchi] = newPatchi++;
+        }
+    }
+    const_cast<pointBoundaryMesh&>(pointBm).reorder(oldToNew, true);
+
+    return pointBm.size() - nPointPatches;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -152,6 +706,12 @@ int main(int argc, char *argv[])
         true  // mark as an advanced option
     );
     argList::addOptionCompat("exclude-patches", {"excludePatches", 2306});
+    argList::addOption
+    (
+        "featureAngle",
+        "angle",
+        "Auto-extract feature edges/points and put into separate point-patches"
+    );
 
     #include "setRootCase.H"
     #include "createTime.H"
@@ -199,9 +759,28 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
+    scalar featureAngle = 180.0;
+    const bool specifiedFeature = args.readIfPresent
+    (
+        "featureAngle",
+        featureAngle
+    );
+
     Info<< "Reading mesh from time " << runTime.value() << endl;
 
     #include "createNamedPolyMesh.H"
+    if (specifiedFeature)
+    {
+        Info<< "Detecting all sharp (>" << featureAngle
+            << " degrees) patch edges." << nl << endl;
+
+        //#include "createNamedPointMesh.H"
+        // Do not read constant/pointMesh - construct from polyMesh only
+        Info<< "Create pointMesh for time = "
+             << runTime.timeName() << Foam::nl << Foam::endl;
+        (void)pointMesh::New(mesh);
+    }
+
 
     // User specified times
     instantList timeDirs = timeSelector::select0(runTime, args);
@@ -275,6 +854,8 @@ int main(int argc, char *argv[])
         // Mesh face and compact zone indx
         DynamicList<label> faceLabels;
         DynamicList<label> compactZones;
+        // Per compact 'zone' index the name and location
+        surfZoneIdentifierList surfZones;
 
         {
             // Collect sizes. Hash on names to handle local-only patches (e.g.
@@ -316,6 +897,17 @@ int main(int argc, char *argv[])
             Pstream::broadcast(compactZoneID);
 
 
+            // Zones
+            surfZones.resize_nocopy(compactZoneID.size());
+            forAllConstIters(compactZoneID, iter)
+            {
+                surfZones[*iter] = surfZoneIdentifier(iter.key(), *iter);
+                Info<< "surfZone " << *iter
+                    <<  " : "      << surfZones[*iter].name()
+                    << endl;
+            }
+
+
             // Rework HashTable into labelList just for speed of conversion
             labelList patchToCompactZone(bMesh.size(), -1);
             labelList faceZoneToCompactZone(bMesh.size(), -1);
@@ -361,7 +953,7 @@ int main(int argc, char *argv[])
 
 
         // Addressing engine for all faces
-        uindirectPrimitivePatch allBoundary
+        const uindirectPrimitivePatch allBoundary
         (
             UIndirectList<face>(mesh.faces(), faceLabels),
             mesh.points()
@@ -399,7 +991,7 @@ int main(int argc, char *argv[])
 
         // Gather all ZoneIDs
         List<labelList> gatheredZones(Pstream::nProcs());
-        gatheredZones[Pstream::myProcNo()].transfer(compactZones);
+        gatheredZones[Pstream::myProcNo()] = compactZones;
         Pstream::gatherList(gatheredZones);
 
         // On master combine all points, faces, zones
@@ -427,16 +1019,6 @@ int main(int argc, char *argv[])
             gatheredZones.clear();
 
 
-            // Zones
-            surfZoneIdentifierList surfZones(compactZoneID.size());
-            forAllConstIters(compactZoneID, iter)
-            {
-                surfZones[*iter] = surfZoneIdentifier(iter.key(), *iter);
-                Info<< "surfZone " << *iter
-                    <<  " : "      << surfZones[*iter].name()
-                    << endl;
-            }
-
             UnsortedMeshedSurface<face> unsortedFace
             (
                 std::move(allPoints),
@@ -463,6 +1045,31 @@ int main(int argc, char *argv[])
 
             sortedFace.write(globalCasePath);
         }
+
+
+        if (specifiedFeature)
+        {
+            // Add edge patches
+            const auto& pMesh = pointMesh::New(mesh);
+
+            const label nAdded = addMeshPointPatches
+            (
+                mesh,
+                pMesh,
+
+                allBoundary,    // all patches together
+                compactZones,   // originating compactZone
+                surfZones,      // per compactZone the index
+
+                featureAngle,
+                featureAngle
+            );
+
+            if (nAdded)
+            {
+                pMesh.boundary().write();
+            }
+        }
     }
 
     Info<< "End\n" << endl;

etc/caseDicts/annotated/topoSetSourcesDict

@@ -440,6 +440,15 @@ pointSet_doc
     }
 
 
+    //- All points of pointpatch
+    {
+        source  patchToPoint;
+        patches ("patch.*");
+        // or
+        patch   somePatch;
+    }
+
+
     //- Copy elements from pointSet
     {
         source  pointToPoint;

src/ODE/ODESolvers/ODESolver/ODESolver.H

@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ODESolver_H
-#define ODESolver_H
+#ifndef Foam_ODESolver_H
+#define Foam_ODESolver_H
 
 #include "ODESystem.H"
 #include "typeInfo.H"
@@ -56,10 +56,9 @@ namespace Foam
 
 class ODESolver
 {
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Reference to ODESystem
         const ODESystem& odes_;
@@ -106,7 +105,7 @@ public:
 
     class stepState
     {
-        public:
+    public:
 
         const bool forward;
         scalar dxTry;
@@ -171,14 +170,14 @@ public:
 
     // Member Functions
 
-        //- Return the number of equations to solve
-        inline label nEqns() const;
+        //- The number of equations to solve
+        label nEqns() const noexcept { return n_; }
 
-        //- Return access to the absolute tolerance field
-        inline scalarField& absTol();
+        //- Access to the absolute tolerance field
+        scalarField& absTol() noexcept { return absTol_; }
 
-        //- Return access to the relative tolerance field
-        inline scalarField& relTol();
+        //- Access to the relative tolerance field
+        scalarField& relTol() noexcept { return relTol_; }
 
         //- Resize the ODE solver
         virtual bool resize() = 0;

src/ODE/ODESolvers/ODESolver/ODESolverI.H

@@ -27,36 +27,19 @@ License
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
-inline Foam::label Foam::ODESolver::nEqns() const
-{
-    return n_;
-}
-
-
-inline Foam::scalarField& Foam::ODESolver::absTol()
-{
-    return absTol_;
-}
-
-
-inline Foam::scalarField& Foam::ODESolver::relTol()
-{
-    return relTol_;
-}
-
-
 template<class Type>
 inline void Foam::ODESolver::resizeField(UList<Type>& f, const label n)
 {
-    f.shallowCopy(UList<Type>(f.begin(), n));
+    // shallowResize
+    f.shallowCopy(f.data(), n);
 }
 
 
 template<class Type>
 inline void Foam::ODESolver::resizeField(UList<Type>& f) const
 {
-    resizeField(f, n_);
+    // shallowResize
+    f.shallowCopy(f.data(), n_);
 }
 
 

src/OpenFOAM/Make/files

@@ -716,8 +716,10 @@ $(pointMeshMapper)/pointPatchMapper.C
 
 pointPatches = $(pointMesh)/pointPatches
 $(pointPatches)/pointPatch/pointPatch.C
+$(pointPatches)/pointPatch/pointPatchNew.C
 $(pointPatches)/facePointPatch/facePointPatch.C
 $(pointPatches)/facePointPatch/facePointPatchNew.C
+$(pointPatches)/meshPointPatch/meshPointPatch.C
 
 basicPointPatches = $(pointPatches)/basic
 $(basicPointPatches)/coupled/coupledPointPatch.C

src/OpenFOAM/containers/PtrLists/PtrListDetail/PtrListDetail.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,7 +30,7 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T>
-Foam::label Foam::Detail::PtrListDetail<T>::count() const noexcept
+Foam::label Foam::Detail::PtrListDetail<T>::count_nonnull() const noexcept
 {
     label n = 0;
 

src/OpenFOAM/containers/PtrLists/PtrListDetail/PtrListDetail.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -95,8 +95,8 @@ public:
         //- Return pointer to element or nullptr for out-of-range access.
         inline T* get(const label i);
 
-        //- Return the number of non-null entries
-        label count() const noexcept;
+        //- The number of non-nullptr entries in the list
+        label count_nonnull() const noexcept;
 
         //- FatalError if any null exists in the list
         inline void checkNonNull() const;

src/OpenFOAM/containers/PtrLists/PtrListDetail/PtrListDetailIO.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,8 +40,8 @@ Foam::Ostream& Foam::Detail::PtrListDetail<T>::write
 {
     const label len = this->size();
 
-    // The net length (after trimming any nullptr)
-    const label netLen = (trimNull ? this->count() : len);
+    // The net length, optionally after trimming any nullptr
+    const label netLen = (trimNull ? this->count_nonnull() : len);
 
     if (!netLen)
     {

src/OpenFOAM/containers/PtrLists/UPtrList/UPtrList.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -257,8 +257,8 @@ public:
         //- Size of the underlying storage.
         inline label capacity() const noexcept;
 
-        //- The number of non-null entries in the list
-        inline label count() const noexcept;
+        //- The number of non-nullptr entries in the list
+        inline label count_nonnull() const noexcept;
 
         //- Reference to the first element of the list
         inline T& front();
@@ -648,6 +648,10 @@ public:
         //- Move append another list to the end of this list.
         //FOAM_DEPRECATED_FOR(2022-10, "push_back()")
         void append(UPtrList<T>&& other) { this->push_back(std::move(other)); }
+
+        //- The number of non-nullptr entries in the list
+        FOAM_DEPRECATED_FOR(2024-01, "count_nonnull()")
+        label count() const noexcept { return count_nonnull(); }
 };
 
 

src/OpenFOAM/containers/PtrLists/UPtrList/UPtrListI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -124,9 +124,9 @@ inline Foam::label Foam::UPtrList<T>::capacity() const noexcept
 
 
 template<class T>
-inline Foam::label Foam::UPtrList<T>::count() const noexcept
+inline Foam::label Foam::UPtrList<T>::count_nonnull() const noexcept
 {
-    return ptrs_.count();
+    return ptrs_.count_nonnull();
 }
 
 

src/OpenFOAM/db/IOstreams/Fstreams/fstreamPointer.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2023 OpenCFD Ltd.
+    Copyright (C) 2020-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -132,7 +132,7 @@ public:
     // Member Functions
 
         //- True if compiled with libz support
-        static bool supports_gz();
+        static bool supports_gz() noexcept;
 
 
     // Access
@@ -276,7 +276,7 @@ public:
     // Member Functions
 
         //- True if compiled with libz support
-        static bool supports_gz();
+        static bool supports_gz() noexcept;
 
 
     // Access

src/OpenFOAM/db/IOstreams/Fstreams/fstreamPointers.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -40,7 +40,7 @@ License
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-bool Foam::ifstreamPointer::supports_gz()
+bool Foam::ifstreamPointer::supports_gz() noexcept
 {
     #ifdef HAVE_LIBZ
     return true;
@@ -50,7 +50,7 @@ bool Foam::ifstreamPointer::supports_gz()
 }
 
 
-bool Foam::ofstreamPointer::supports_gz()
+bool Foam::ofstreamPointer::supports_gz() noexcept
 {
     #ifdef HAVE_LIBZ
     return true;
@@ -71,7 +71,7 @@ Foam::ifstreamPointer::ifstreamPointer
     IOstreamOption streamOpt  // Currently unused
 )
 :
-    ptr_(nullptr)
+    ptr_()
 {
     open(pathname, streamOpt);
 }
@@ -82,7 +82,7 @@ Foam::ifstreamPointer::ifstreamPointer
     const fileName& pathname
 )
 :
-    ptr_(nullptr)
+    ptr_()
 {
     open(pathname);
 }
@@ -110,7 +110,7 @@ Foam::ofstreamPointer::ofstreamPointer
     const bool atomic
 )
 :
-    ptr_(nullptr),
+    ptr_(),
     atomic_(atomic)
 {
     std::ios_base::openmode mode

src/OpenFOAM/db/IOstreams/IOstreams/IOstream.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -468,6 +468,41 @@ inline IOstream& scientific(IOstream& io)
 namespace Detail
 {
 
+/*---------------------------------------------------------------------------*\
+                   Class Detail::StreamAllocator Declaration
+\*---------------------------------------------------------------------------*/
+
+//- A wrapper to hold a std::stream type for OpenFOAM wrapped streams.
+//- This is necessary since the OpenFOAM streams hold a reference to
+//- the normal std::stream
+template<class StreamType>
+class StreamAllocator
+{
+protected:
+
+    // Protected Data
+
+        //- The std::stream
+        StreamType stream_;
+
+
+    // Constructors
+
+        //- Default construct (empty)
+        StreamAllocator() = default;
+};
+
+
+} // End namespace Detail
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Functions/Algorithms
+
+namespace Detail
+{
+
 //- Termination for input looping (no-op)
 template<class IS> inline void inputLoop(IS&) {}
 

src/OpenFOAM/db/IOstreams/StringStreams/StringStream.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,60 +46,6 @@ SourceFiles
 namespace Foam
 {
 
-namespace Detail
-{
-
-/*---------------------------------------------------------------------------*\
-                Class Detail::StringStreamAllocator Declaration
-\*---------------------------------------------------------------------------*/
-
-//- Allocator for variants of a std stringstream
-template<class StreamType>
-class StringStreamAllocator
-{
-protected:
-
-    // Protected Member Data
-
-        //- The stream type
-        typedef StreamType stream_type;
-
-        //- The input/output stream.
-        stream_type stream_;
-
-
-    // Constructors
-
-        //- Default construct
-        StringStreamAllocator() = default;
-
-        //- Copy construct from string
-        StringStreamAllocator(const std::string& s)
-        :
-            stream_(s)
-        {}
-
-
-public:
-
-    // Member Functions
-
-        //- Get the string - as Foam::string rather than std::string
-        Foam::string str() const
-        {
-            return Foam::string(stream_.str());
-        }
-
-        //- Set the string
-        void str(const std::string& s)
-        {
-            stream_.str(s);
-        }
-};
-
-} // End namespace Detail
-
-
 /*---------------------------------------------------------------------------*\
                         Class IStringStream Declaration
 \*---------------------------------------------------------------------------*/
@@ -107,10 +53,12 @@ public:
 //- Input from string buffer, using a ISstream. Always UNCOMPRESSED.
 class IStringStream
 :
-    public Detail::StringStreamAllocator<std::istringstream>,
-    public ISstream
+    public Foam::Detail::StreamAllocator<std::istringstream>,
+    public Foam::ISstream
 {
-    typedef Detail::StringStreamAllocator<std::istringstream> allocator_type;
+    typedef
+        Foam::Detail::StreamAllocator<std::istringstream>
+        allocator_type;
 
 public:
 
@@ -133,9 +81,11 @@ public:
             IOstreamOption streamOpt = IOstreamOption()
         )
         :
-            allocator_type(s),
+            allocator_type(),
             ISstream(stream_, "input", streamOpt.format(), streamOpt.version())
-        {}
+        {
+            stream_.str(s);
+        }
 
         //- Construct from char*
         explicit IStringStream
@@ -144,20 +94,32 @@ public:
             IOstreamOption streamOpt = IOstreamOption()
         )
         :
-            allocator_type(s),
+            allocator_type(),
             ISstream(stream_, "input", streamOpt.format(), streamOpt.version())
-        {}
+        {
+            stream_.str(s);
+        }
 
         //- Copy construct, copies content and format
         IStringStream(const IStringStream& str)
         :
-            allocator_type(str.str()),
+            allocator_type(),
             ISstream(stream_, str.name(), static_cast<IOstreamOption>(str))
-        {}
+        {
+            stream_.str(str.str());
+        }
 
 
     // Member Functions
 
+        //- Get the string.
+        //- As Foam::string instead of std::string (may change in future)
+        Foam::string str() const { return Foam::string(stream_.str()); }
+
+        //- Set the string
+        void str(const std::string& s) { stream_.str(s); }
+
+
         //- Reset the input buffer and rewind the stream
         virtual void reset(const std::string& s)
         {
@@ -220,10 +182,12 @@ public:
 //- Output to string buffer, using a OSstream. Always UNCOMPRESSED.
 class OStringStream
 :
-    public Detail::StringStreamAllocator<std::ostringstream>,
-    public OSstream
+    public Foam::Detail::StreamAllocator<std::ostringstream>,
+    public Foam::OSstream
 {
-    typedef Detail::StringStreamAllocator<std::ostringstream> allocator_type;
+    typedef
+        Foam::Detail::StreamAllocator<std::ostringstream>
+        allocator_type;
 
 public:
 
@@ -242,13 +206,23 @@ public:
         //- Copy construct, copies content and format
         OStringStream(const OStringStream& str)
         :
-            allocator_type(str.str()),
+            allocator_type(),
             OSstream(stream_, str.name(), static_cast<IOstreamOption>(str))
-        {}
+        {
+            stream_.str(str.str());
+        }
 
 
     // Member Functions
 
+        //- Get the string.
+        //- As Foam::string instead of std::string (may change in future)
+        Foam::string str() const { return Foam::string(stream_.str()); }
+
+        //- Set the string
+        void str(const std::string& s) { stream_.str(s); }
+
+
         //- Reset the output buffer and rewind the stream
         void reset()
         {

src/OpenFOAM/db/IOstreams/Tstreams/ITstream.C

@@ -91,13 +91,14 @@ Foam::ITstream& Foam::ITstream::empty_stream()
 }
 
 
-Foam::tokenList Foam::ITstream::parse
+Foam::tokenList Foam::ITstream::parse_chars
 (
-    const UList<char>& input,
+    const char* s,
+    size_t nbytes,
     IOstreamOption streamOpt
 )
 {
-    ISpanStream is(input, streamOpt);
+    ISpanStream is(s, nbytes, streamOpt);
 
     tokenList tokens;
     parseStream(is, tokens);
@@ -105,31 +106,14 @@ Foam::tokenList Foam::ITstream::parse
 }
 
 
-Foam::tokenList Foam::ITstream::parse
-(
-    const std::string& input,
-    IOstreamOption streamOpt
-)
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::ITstream::reset(const char* input, size_t nbytes)
 {
-    ISpanStream is(input, streamOpt);
+    ISpanStream is(input, nbytes, static_cast<IOstreamOption>(*this));
 
-    tokenList tokens;
-    parseStream(is, tokens);
-    return tokens;
-}
-
-
-Foam::tokenList Foam::ITstream::parse
-(
-    const char* input,
-    IOstreamOption streamOpt
-)
-{
-    ISpanStream is(input, strlen(input), streamOpt);
-
-    tokenList tokens;
-    parseStream(is, tokens);
-    return tokens;
+    parseStream(is, static_cast<tokenList&>(*this));
+    ITstream::seek(0);  // rewind(), but bypasss virtual
 }
 
 
@@ -254,10 +238,7 @@ Foam::ITstream::ITstream
 :
     ITstream(streamOpt, name)
 {
-    ISpanStream is(input, streamOpt);
-
-    parseStream(is, static_cast<tokenList&>(*this));
-    ITstream::seek(0);  // rewind(), but bypasss virtual
+    reset(input.cdata(), input.size_bytes());
 }
 
 
@@ -270,10 +251,7 @@ Foam::ITstream::ITstream
 :
     ITstream(streamOpt, name)
 {
-    ISpanStream is(input, streamOpt);
-
-    parseStream(is, static_cast<tokenList&>(*this));
-    ITstream::seek(0);  // rewind(), but bypasss virtual
+    reset(input.data(), input.size());
 }
 
 
@@ -286,10 +264,7 @@ Foam::ITstream::ITstream
 :
     ITstream(streamOpt, name)
 {
-    ISpanStream is(input, strlen(input), streamOpt);
-
-    parseStream(is, static_cast<tokenList&>(*this));
-    ITstream::seek(0);  // rewind(), but bypasss virtual
+    reset(input, strlen(input));
 }
 
 

src/OpenFOAM/db/IOstreams/Tstreams/ITstream.H

@@ -77,6 +77,18 @@ class ITstream
         //  but leave any excess capacity (ie, like reserve).
         void reserveCapacity(const label newCapacity);
 
+        //- Convert input sequence into a list of tokens,
+        static tokenList parse_chars
+        (
+            const char* s,
+            size_t nbytes,
+            IOstreamOption streamOpt
+        );
+
+        //- Convert input sequence into a list of tokens,
+        //- using the existing stream format. Rewinds the stream
+        void reset(const char* input, size_t nbytes);
+
         //- Failsafe read-access to token at specified location
         //- or undefinedToken
         inline const token& peekNoFail(const label i) const
@@ -158,6 +170,47 @@ public:
             const string& name = "input"
         );
 
+        #if __cplusplus >= 201703L
+        //- Construct token list by parsing the input character sequence
+        //  Uses static parse function internally.
+        explicit ITstream
+        (
+            std::string_view s,
+            IOstreamOption streamOpt = IOstreamOption()
+        )
+        :
+            ITstream(streamOpt)
+        {
+            reset(s.data(), s.size());
+        }
+        #endif
+
+        //- Construct token list by parsing the input character sequence
+        //  Uses static parse function internally.
+        explicit ITstream
+        (
+            stdFoam::span<char> s,
+            IOstreamOption streamOpt = IOstreamOption()
+        )
+        :
+            ITstream(streamOpt)
+        {
+            reset(s.data(), s.size());
+        }
+
+        //- Construct token list by parsing the input character sequence
+        //  Uses static parse function internally.
+        explicit ITstream
+        (
+            stdFoam::span<const char> s,
+            IOstreamOption streamOpt = IOstreamOption()
+        )
+        :
+            ITstream(streamOpt)
+        {
+            reset(s.data(), s.size());
+        }
+
 
     // Additional constructors
 
@@ -207,7 +260,10 @@ public:
         (
             const UList<char>& input,
             IOstreamOption streamOpt = IOstreamOption()
-        );
+        )
+        {
+            return parse_chars(input.cdata(), input.size(), streamOpt);
+        }
 
         //- Create token list by parsing the input string
         //- until no good tokens remain.
@@ -215,7 +271,10 @@ public:
         (
             const std::string& input,
             IOstreamOption streamOpt = IOstreamOption()
-        );
+        )
+        {
+            return parse_chars(input.data(), input.size(), streamOpt);
+        }
 
         //- Create token list by parsing the input character sequence
         //- until no good tokens remain.
@@ -223,7 +282,45 @@ public:
         (
             const char* input,
             IOstreamOption streamOpt = IOstreamOption()
-        );
+        )
+        {
+            return parse_chars(input, strlen(input), streamOpt);
+        }
+
+        #if __cplusplus >= 201703L
+        //- Create token list by parsing the input character sequence
+        //- until no good tokens remain.
+        static tokenList parse
+        (
+            std::string_view s,
+            IOstreamOption streamOpt = IOstreamOption()
+        )
+        {
+            return parse_chars(s.data(), s.size(), streamOpt);
+        }
+        #endif
+
+        //- Create token list by parsing the input character sequence
+        //- until no good tokens remain.
+        static tokenList parse
+        (
+            stdFoam::span<char> s,
+            IOstreamOption streamOpt = IOstreamOption()
+        )
+        {
+            return parse_chars(s.data(), s.size(), streamOpt);
+        }
+
+        //- Create token list by parsing the input character sequence
+        //- until no good tokens remain.
+        static tokenList parse
+        (
+            stdFoam::span<const char> s,
+            IOstreamOption streamOpt = IOstreamOption()
+        )
+        {
+            return parse_chars(s.data(), s.size(), streamOpt);
+        }
 
 
     // Member Functions

src/OpenFOAM/db/IOstreams/hashes/OSHA1stream.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2019-2023 OpenCFD Ltd.
+    Copyright (C) 2019-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,8 +35,8 @@ Description
 #ifndef Foam_OSHA1stream_H
 #define Foam_OSHA1stream_H
 
-#include "OSstream.H"
 #include "SHA1.H"
+#include "OSstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -114,33 +114,6 @@ public:
 };
 
 
-namespace Detail
-{
-
-/*---------------------------------------------------------------------------*\
-                Class Detail::OSHA1streamAllocator Declaration
-\*---------------------------------------------------------------------------*/
-
-//- An allocator for holding Foam::osha1stream
-class OSHA1streamAllocator
-{
-protected:
-
-    // Protected Data
-
-        //- The output stream
-        Foam::osha1stream stream_;
-
-
-    // Constructors
-
-        //- Default construct
-        OSHA1streamAllocator() = default;
-};
-
-} // End namespace Detail
-
-
 /*---------------------------------------------------------------------------*\
                          Class OSHA1stream Declaration
 \*---------------------------------------------------------------------------*/
@@ -148,10 +121,12 @@ protected:
 //- The output stream for calculating SHA1 digests
 class OSHA1stream
 :
-    public Detail::OSHA1streamAllocator,
-    public OSstream
+    public Foam::Detail::StreamAllocator<Foam::osha1stream>,
+    public Foam::OSstream
 {
-    typedef Detail::OSHA1streamAllocator allocator_type;
+    typedef
+        Foam::Detail::StreamAllocator<Foam::osha1stream>
+        allocator_type;
 
 public:
 

src/OpenFOAM/db/IOstreams/memory/ICharStream.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2023 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,8 +42,6 @@ See Also
 #define Foam_ICharStream_H
 
 #include "ISpanStream.H"
-#include "List.H"
-#include "DynamicList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -205,33 +203,6 @@ public:
 };
 
 
-namespace Detail
-{
-
-/*---------------------------------------------------------------------------*\
-                Class Detail::ICharStreamAllocator Declaration
-\*---------------------------------------------------------------------------*/
-
-//- An allocator for holding Foam::icharstream
-class ICharStreamAllocator
-{
-protected:
-
-    // Protected Data
-
-        //- The stream
-        Foam::icharstream stream_;
-
-
-    // Constructors
-
-        //- Default construct
-        ICharStreamAllocator() = default;
-};
-
-} // End namespace Detail
-
-
 /*---------------------------------------------------------------------------*\
                          Class ICharStream Declaration
 \*---------------------------------------------------------------------------*/
@@ -239,10 +210,12 @@ protected:
 //- An ISstream with internal List storage. Always UNCOMPRESSED.
 class ICharStream
 :
-    public Detail::ICharStreamAllocator,
+    public Foam::Detail::StreamAllocator<Foam::icharstream>,
     public Foam::ISstream
 {
-    typedef Detail::ICharStreamAllocator allocator_type;
+    typedef
+        Foam::Detail::StreamAllocator<Foam::icharstream>
+        allocator_type;
 
 public:
 

src/OpenFOAM/db/IOstreams/memory/ISpanStream.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2023 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,7 +65,6 @@ See Also
 #define Foam_ISpanStream_H
 
 #include "memoryStreamBuffer.H"
-#include "UList.H"
 #include "ISstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -213,49 +212,23 @@ public:
         //- Some information about the input buffer position/capacity
         void debug_info(Ostream& os) const
         {
-            os  << "get="
-                << input_pos() << '/' << capacity();
+            os  << "get=" << input_pos() << '/' << capacity();
         }
 };
 
 
-namespace Detail
-{
-
-/*---------------------------------------------------------------------------*\
-                Class Detail::ISpanStreamAllocator Declaration
-\*---------------------------------------------------------------------------*/
-
-//- An allocator for holding Foam::ispanstream
-class ISpanStreamAllocator
-{
-protected:
-
-    // Protected Data
-
-        //- The stream
-        Foam::ispanstream stream_;
-
-
-    // Constructors
-
-        //- Default construct (empty)
-        ISpanStreamAllocator() = default;
-};
-
-} // End namespace Detail
-
-
 /*---------------------------------------------------------------------------*\
                          Class ISpanStream Declaration
 \*---------------------------------------------------------------------------*/
 
 class ISpanStream
 :
-    public Detail::ISpanStreamAllocator,
+    public Foam::Detail::StreamAllocator<Foam::ispanstream>,
     public Foam::ISstream
 {
-    typedef Detail::ISpanStreamAllocator allocator_type;
+    typedef
+        Foam::Detail::StreamAllocator<Foam::ispanstream>
+        allocator_type;
 
 public:
 

src/OpenFOAM/db/IOstreams/memory/OCharStream.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2017-2023 OpenCFD Ltd.
+    Copyright (C) 2017-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -109,7 +109,13 @@ public:
         //- The current output position within the buffer (tellp)
         std::streampos output_pos() const
         {
-            return (buffer_type::span_tellp());
+            return buffer_type::span_tellp();
+        }
+
+        //- The number of bytes outputted
+        std::streamsize count() const
+        {
+            return buffer_type::size_bytes();
         }
 
         //- The put buffer capacity
@@ -189,33 +195,6 @@ public:
 };
 
 
-namespace Detail
-{
-
-/*---------------------------------------------------------------------------*\
-                Class Detail::OCharStreamAllocator Declaration
-\*---------------------------------------------------------------------------*/
-
-//- An allocator for holding Foam::ocharstream
-class OCharStreamAllocator
-{
-protected:
-
-    // Protected Data
-
-        //- The stream
-        Foam::ocharstream stream_;
-
-
-    // Constructors
-
-        //- Default construct - empty
-        OCharStreamAllocator() = default;
-};
-
-} // End namespace Detail
-
-
 /*---------------------------------------------------------------------------*\
                          Class OCharStream Declaration
 \*---------------------------------------------------------------------------*/
@@ -223,10 +202,12 @@ protected:
 //- An OSstream with internal List storage
 class OCharStream
 :
-    public Detail::OCharStreamAllocator,
+    public Foam::Detail::StreamAllocator<Foam::ocharstream>,
     public Foam::OSstream
 {
-    typedef Detail::OCharStreamAllocator allocator_type;
+    typedef
+        Foam::Detail::StreamAllocator<Foam::ocharstream>
+        allocator_type;
 
 public:
 
@@ -277,8 +258,11 @@ public:
         //- The current output position within the buffer (tellp)
         std::streampos output_pos() const { return stream_.output_pos(); }
 
-        //- The current output size. Same as tellp(), output_pos()
-        label size() const { return label(stream_.output_pos()); }
+        //- The number of bytes outputted
+        std::streamsize count() const { return stream_.count(); }
+
+        //- The current output size. Same as count(), output_pos(), tellp().
+        label size() const { return label(stream_.count()); }
 
         //- The put buffer capacity
         std::streamsize capacity() const { return stream_.capacity(); }

src/OpenFOAM/db/IOstreams/memory/OCountStream.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2023 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,7 +35,6 @@ Description
 #define Foam_OScountStream_H
 
 #include "OSstream.H"
-#include <iostream>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -126,7 +125,7 @@ class ocountstream
 
         // Member Functions
 
-            //- The number of bytes counted.
+            //- The number of bytes counted
             std::streamsize count() const noexcept { return size_; }
 
             //- Reset the count
@@ -155,7 +154,7 @@ public:
         //- This hides both signatures of std::basic_ios::rdbuf()
         countbuf* rdbuf() { return &buf_; }
 
-        //- \return The number of bytes counted
+        //- The number of bytes counted
         std::streamsize count() const noexcept { return buf_.count(); }
 
         //- Reset the count
@@ -173,33 +172,6 @@ public:
 };
 
 
-namespace Detail
-{
-
-/*---------------------------------------------------------------------------*\
-                Class Detail::OCountStreamAllocator Declaration
-\*---------------------------------------------------------------------------*/
-
-//- An allocator for holding Foam::ocountstream
-class OCountStreamAllocator
-{
-protected:
-
-    // Protected Data
-
-        //- The output stream
-        Foam::ocountstream stream_;
-
-
-    // Constructors
-
-        //- Default construct
-        OCountStreamAllocator() = default;
-};
-
-} // End namespace Detail
-
-
 /*---------------------------------------------------------------------------*\
                         Class OCountStream Declaration
 \*---------------------------------------------------------------------------*/
@@ -207,10 +179,12 @@ protected:
 //- An output stream for calculating byte counts
 class OCountStream
 :
-    public Detail::OCountStreamAllocator,
-    public OSstream
+    public Foam::Detail::StreamAllocator<Foam::ocountstream>,
+    public Foam::OSstream
 {
-    typedef Detail::OCountStreamAllocator allocator_type;
+    typedef
+        Foam::Detail::StreamAllocator<Foam::ocountstream>
+        allocator_type;
 
 public:
 
@@ -238,10 +212,10 @@ public:
 
     // Member Functions
 
-        //- \return The number of bytes counted
+        //- The number of bytes counted
         std::streamsize count() const noexcept { return stream_.count(); }
 
-        //- \return The number of bytes counted
+        //- The number of bytes counted
         std::streamsize size() const noexcept { return stream_.count(); }
 
         //- Reset the count

src/OpenFOAM/db/IOstreams/memory/OSpanStream.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2023 OpenCFD Ltd.
+    Copyright (C) 2016-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,7 +83,6 @@ See Also
 #define Foam_OSpanStream_H
 
 #include "memoryStreamBuffer.H"
-#include "DynamicList.H"
 #include "OSstream.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -141,6 +140,12 @@ public:
             return buffer_type::span_tellp();
         }
 
+        //- The number of bytes outputted
+        std::streamsize count() const
+        {
+            return buffer_type::size_bytes();
+        }
+
         //- The put buffer capacity
         std::streamsize capacity() const
         {
@@ -199,49 +204,23 @@ public:
         //- Some information about the output buffer position/capacity
         void debug_info(Ostream& os) const
         {
-            os  << "put="
-                << output_pos() << '/' << capacity();
+            os  << "put=" << output_pos() << '/' << capacity();
         }
 };
 
 
-namespace Detail
-{
-
-/*---------------------------------------------------------------------------*\
-                Class Detail::OSpanStreamAllocator Declaration
-\*---------------------------------------------------------------------------*/
-
-//- An allocator for holding Foam::ospanstream
-class OSpanStreamAllocator
-{
-protected:
-
-    // Protected Data
-
-        //- The stream
-        Foam::ospanstream stream_;
-
-
-    // Constructors
-
-        //- Default construct (empty)
-        OSpanStreamAllocator() = default;
-};
-
-} // End namespace Detail
-
-
 /*---------------------------------------------------------------------------*\
                          Class OSpanStream Declaration
 \*---------------------------------------------------------------------------*/
 
 class OSpanStream
 :
-    public Detail::OSpanStreamAllocator,
+    public Foam::Detail::StreamAllocator<Foam::ospanstream>,
     public Foam::OSstream
 {
-    typedef Detail::OSpanStreamAllocator allocator_type;
+    typedef
+        Foam::Detail::StreamAllocator<Foam::ospanstream>
+        allocator_type;
 
 public:
 
@@ -314,8 +293,11 @@ public:
         //- The current output position within the buffer (tellp)
         std::streampos output_pos() const { return stream_.output_pos(); }
 
-        //- The current output size. Same as tellp(), output_pos()
-        label size() const { return label(stream_.output_pos()); }
+        //- The number of bytes outputted
+        std::streamsize count() const { return stream_.count(); }
+
+        //- The current output size. Same as count(), output_pos(), tellp().
+        label size() const { return label(stream_.count()); }
 
         //- The put buffer capacity
         std::streamsize capacity() const { return stream_.capacity(); }

src/OpenFOAM/include/createNamedPointMesh.H

@@ -0,0 +1,10 @@
+Foam::Info
+    << "Create pointMesh for time = "
+    << runTime.timeName() << Foam::nl << Foam::endl;
+
+// Register pointMesh on the database
+const Foam::pointMesh& pMesh = pointMesh::New
+(
+    mesh,
+    Foam::IOobject::READ_IF_PRESENT
+);

src/OpenFOAM/matrices/LUscalarMatrix/procLduMatrix.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,22 +43,11 @@ Foam::procLduMatrix::procLduMatrix
     lowerAddr_(ldum.lduAddr().lowerAddr()),
     diag_(ldum.diag()),
     upper_(ldum.upper()),
-    lower_(ldum.lower())
+    lower_(ldum.lower()),
+    interfaces_(interfaces.count_nonnull())
 {
     label nInterfaces = 0;
 
-    forAll(interfaces, i)
-    {
-        if (interfaces.set(i))
-        {
-            nInterfaces++;
-        }
-    }
-
-    interfaces_.setSize(nInterfaces);
-
-    nInterfaces = 0;
-
     forAll(interfaces, i)
     {
         if (interfaces.set(i))
@@ -73,7 +63,6 @@ Foam::procLduMatrix::procLduMatrix
             );
         }
     }
-
 }
 
 

src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,9 +33,74 @@ License
 #include "PstreamBuffers.H"
 #include "lduSchedule.H"
 #include "globalMeshData.H"
+#include "processorPointPatch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(pointBoundaryMesh, 0);
+}
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+bool Foam::pointBoundaryMesh::hasGroupIDs() const
+{
+    if (groupIDsPtr_)
+    {
+        // Use existing cache
+        return !groupIDsPtr_->empty();
+    }
+
+    const auto& patches = *this;
+
+    for (const auto& p : patches)
+    {
+        if (!p.inGroups().empty())
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+void Foam::pointBoundaryMesh::calcGroupIDs() const
+{
+    if (groupIDsPtr_)
+    {
+        return;  // Or FatalError
+    }
+
+    groupIDsPtr_.emplace(16);
+    auto& groupLookup = *groupIDsPtr_;
+
+    const auto& patches = *this;
+
+    forAll(patches, patchi)
+    {
+        for (const word& groupName : patches[patchi].inGroups())
+        {
+            groupLookup(groupName).push_back(patchi);
+        }
+    }
+
+    // Remove groups that clash with patch names
+    forAll(patches, patchi)
+    {
+        if (groupLookup.erase(patches[patchi].name()))
+        {
+            WarningInFunction
+                << "Removed group '" << patches[patchi].name()
+                << "' which clashes with patch " << patchi
+                << " of the same name."
+                << endl;
+        }
+    }
+}
+
+
 void Foam::pointBoundaryMesh::addPatches(const polyBoundaryMesh& pbm)
 {
     // Set boundary patches
@@ -43,7 +108,7 @@ void Foam::pointBoundaryMesh::addPatches(const polyBoundaryMesh& pbm)
 
     patches.resize_null(pbm.size());
 
-    forAll(patches, patchi)
+    forAll(pbm, patchi)
     {
         // NB: needs ptr() to get *pointPatch instead of *facePointPatch
         patches.set(patchi, facePointPatch::New(pbm[patchi], *this).ptr());
@@ -60,21 +125,301 @@ Foam::pointBoundaryMesh::pointBoundaryMesh
 )
 :
     pointPatchList(),
+    regIOobject
+    (
+        IOobject
+        (
+            "boundary",
+            //m.thisDb().time().findInstance(m.meshDir(), "boundary"),
+            pbm.mesh().facesInstance(),
+            polyMesh::meshSubDir/pointMesh::meshSubDir,
+            m.thisDb(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false                   // Avoid conflict with polyMesh/boundary
+        )
+    ),
     mesh_(m)
 {
     addPatches(pbm);
+
+    if (debug)
+    {
+        pointPatchList& Patches = *this;
+
+        Pout<< "pointBoundaryMesh::pointBoundaryMesh"
+            << "(const pointMesh&, const polyBoundaryMesh&): "
+            << "constructed pointBoundaryMesh:" << endl;
+        Pout<< incrIndent;
+        for (const auto& pp : Patches)
+        {
+            Pout<< indent
+                << "index:" << pp.index() << " patch:" << pp.name()
+                << " type:" << pp.type() << endl;
+        }
+        Pout<< decrIndent;
+    }
+}
+
+
+Foam::pointBoundaryMesh::pointBoundaryMesh
+(
+    const IOobject& io,
+    const pointMesh& m,
+    const polyBoundaryMesh& pbm
+)
+:
+    pointPatchList(),
+    regIOobject
+    (
+        IOobject
+        (
+            "boundary",
+            io.instance(),
+            polyMesh::meshSubDir/pointMesh::meshSubDir,
+            io.db(),
+            io.readOpt(),
+            io.writeOpt(),
+            false   //io.registerObject()     // or always set to false?
+        )
+    ),
+    mesh_(m)
+{
+    pointPatchList& Patches = *this;
+
+    if (isReadRequired() || (isReadOptional() && headerOk()))
+    {
+        if (readOpt() == IOobject::MUST_READ_IF_MODIFIED)
+        {
+            WarningInFunction
+                << "Specified IOobject::MUST_READ_IF_MODIFIED but class"
+                << " does not support automatic rereading."
+                << endl;
+        }
+
+        if (debug)
+        {
+            Pout<< "pointBoundaryMesh::pointBoundaryMesh"
+                << "(const IOobject&, const pointMesh&,"
+                << " const polyBoundaryMesh&): "
+                << "Constructing from boundary file " << objectRelPath()
+                << endl;
+        }
+
+        // Read pointPatchList
+        Istream& is = readStream(typeName);
+
+        PtrList<entry> patchEntries(is);
+        Patches.setSize(patchEntries.size());
+
+        forAll(Patches, patchi)
+        {
+            // Try construct-from-dictionary first
+            const word& name = patchEntries[patchi].keyword();
+
+            autoPtr<pointPatch> pPtr
+            (
+                pointPatch::New
+                (
+                    name,
+                    patchEntries[patchi].dict(),
+                    patchi,
+                    *this
+                )
+            );
+
+            if (!pPtr)
+            {
+                const label polyPatchi = pbm.findPatchID(name, false);
+                // Try as facePointPatch from polyPatch
+                pPtr = facePointPatch::New(pbm[polyPatchi], *this);
+                pPtr->index() = patchi;
+            }
+
+            Patches.set(patchi, pPtr);
+        }
+
+        // Check state of IOstream
+        is.check
+        (
+            "pointBoundaryMesh::pointBoundaryMesh"
+            "(const IOobject&, const pointMesh&,"
+            " const polyBoundaryMesh&)"
+        );
+
+        close();
+    }
+    else
+    {
+        if (debug)
+        {
+            Pout<< "pointBoundaryMesh::pointBoundaryMesh"
+                << "(const IOobject&, const pointMesh&,"
+                << " const polyBoundaryMesh&): "
+                << "Constructing from polyBoundaryMesh only"
+                << endl;
+        }
+
+        addPatches(pbm);
+    }
+
+
+    if (debug)
+    {
+        Pout<< "pointBoundaryMesh::pointBoundaryMesh"
+            << "(const IOobject&, const pointMesh&, const polyBoundaryMesh&): "
+            << "constructed pointBoundaryMesh:" << endl;
+        Pout<< incrIndent;
+        for (const auto& pp : Patches)
+        {
+            Pout<< indent
+                << "index:" << pp.index() << " patch:" << pp.name()
+                << " type:" << pp.type() << endl;
+        }
+        Pout<< decrIndent;
+    }
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::label Foam::pointBoundaryMesh::nNonProcessor() const
+{
+    const pointPatchList& patches = *this;
+
+    label count = 0;
+
+    for (const auto& p : patches)
+    {
+        if (isA<processorPointPatch>(p))
+        {
+            break;
+        }
+
+        ++count;
+    }
+
+    return count;
+}
+
+
+Foam::label Foam::pointBoundaryMesh::nProcessorPatches() const
+{
+    const pointPatchList& patches = *this;
+
+    label count = 0;
+
+    for (const auto& p : patches)
+    {
+        if (isA<processorPointPatch>(p))
+        {
+            ++count;
+        }
+    }
+
+    return count;
+}
+
+
+Foam::wordList Foam::pointBoundaryMesh::names() const
+{
+    return PtrListOps::get<word>(*this, nameOp<pointPatch>());
+}
+
+
+Foam::wordList Foam::pointBoundaryMesh::types() const
+{
+    return PtrListOps::get<word>(*this, typeOp<pointPatch>());
+}
+
+
+Foam::wordList Foam::pointBoundaryMesh::physicalTypes() const
+{
+    return
+        PtrListOps::get<word>
+        (
+            *this,
+            [](const pointPatch& p) { return p.physicalType(); }
+        );
+}
+
+
+const Foam::HashTable<Foam::labelList>&
+Foam::pointBoundaryMesh::groupPatchIDs() const
+{
+    if (!groupIDsPtr_)
+    {
+        calcGroupIDs();
+    }
+
+    return *groupIDsPtr_;
+}
+
+
 Foam::labelList Foam::pointBoundaryMesh::indices
 (
     const wordRe& matcher,
     const bool useGroups
 ) const
 {
-    return mesh().boundaryMesh().indices(matcher, useGroups);
+    if (matcher.empty())
+    {
+        return labelList();
+    }
+
+    // Only check groups if requested and they exist
+    const bool checkGroups = (useGroups && this->hasGroupIDs());
+
+    labelHashSet ids(0);
+
+    if (matcher.isPattern())
+    {
+        if (checkGroups)
+        {
+            const auto& groupLookup = groupPatchIDs();
+            forAllConstIters(groupLookup, iter)
+            {
+                if (matcher(iter.key()))
+                {
+                    // Add patch ids associated with the group
+                    ids.insert(iter.val());
+                }
+            }
+        }
+
+        if (ids.empty())
+        {
+            return PtrListOps::findMatching(*this, matcher);
+        }
+        else
+        {
+            ids.insert(PtrListOps::findMatching(*this, matcher));
+        }
+    }
+    else
+    {
+        // Literal string.
+        // Special version of above for reduced memory footprint.
+
+        const label patchId = PtrListOps::firstMatching(*this, matcher);
+
+        if (patchId >= 0)
+        {
+            return labelList(one{}, patchId);
+        }
+        else if (checkGroups)
+        {
+            const auto iter = groupPatchIDs().cfind(matcher);
+
+            if (iter.good())
+            {
+                // Hash ids associated with the group
+                ids.insert(iter.val());
+            }
+        }
+    }
+
+    return ids.sortedToc();
 }
 
 
@@ -84,7 +429,43 @@ Foam::labelList Foam::pointBoundaryMesh::indices
     const bool useGroups
 ) const
 {
-    return mesh().boundaryMesh().indices(matcher, useGroups);
+    if (matcher.empty())
+    {
+        return labelList();
+    }
+    else if (matcher.size() == 1)
+    {
+        return this->indices(matcher.front(), useGroups);
+    }
+
+    labelHashSet ids(0);
+
+    // Only check groups if requested and they exist
+    if (useGroups && this->hasGroupIDs())
+    {
+        ids.reserve(this->size());
+
+        const auto& groupLookup = groupPatchIDs();
+        forAllConstIters(groupLookup, iter)
+        {
+            if (matcher(iter.key()))
+            {
+                // Add patch ids associated with the group
+                ids.insert(iter.val());
+            }
+        }
+    }
+
+    if (ids.empty())
+    {
+        return PtrListOps::findMatching(*this, matcher);
+    }
+    else
+    {
+        ids.insert(PtrListOps::findMatching(*this, matcher));
+    }
+
+    return ids.sortedToc();
 }
 
 
@@ -95,13 +476,91 @@ Foam::labelList Foam::pointBoundaryMesh::indices
     const bool useGroups
 ) const
 {
-    return mesh().boundaryMesh().indices(select, ignore, useGroups);
+    //return mesh().boundaryMesh().indices(select, ignore, useGroups);
+    if (ignore.empty())
+    {
+        return this->indices(select, useGroups);
+    }
+
+    const wordRes::filter matcher(select, ignore);
+
+    labelHashSet ids(0);
+
+    // Only check groups if requested and they exist
+    if (useGroups && this->hasGroupIDs())
+    {
+        ids.reserve(this->size());
+
+        const auto& groupLookup = groupPatchIDs();
+        forAllConstIters(groupLookup, iter)
+        {
+            if (matcher(iter.key()))
+            {
+                // Add patch ids associated with the group
+                ids.insert(iter.val());
+            }
+        }
+    }
+
+    if (ids.empty())
+    {
+        return PtrListOps::findMatching(*this, matcher);
+    }
+    else
+    {
+        ids.insert(PtrListOps::findMatching(*this, matcher));
+    }
+
+    return ids.sortedToc();
 }
 
 
-Foam::label Foam::pointBoundaryMesh::findPatchID(const word& patchName) const
+Foam::label Foam::pointBoundaryMesh::findPatchID
+(
+    const word& patchName,
+    bool allowNotFound
+) const
 {
-    return mesh().boundaryMesh().findPatchID(patchName);
+    //return mesh().boundaryMesh().findPatchID(patchName);
+    if (patchName.empty())
+    {
+        return -1;
+    }
+
+    const label patchId = PtrListOps::firstMatching(*this, patchName);
+
+    if (patchId >= 0)
+    {
+        return patchId;
+    }
+
+    if (!allowNotFound)
+    {
+        FatalErrorInFunction
+            << "Patch '" << patchName << "' not found. "
+            << "Available patch names";
+
+        if (polyMesh::defaultRegion != mesh_.name())
+        {
+            FatalError
+                << " in region '" << mesh_.name() << "'";
+        }
+
+        FatalError
+            << " include: " << names() << endl
+            << exit(FatalError);
+    }
+
+    // Patch not found
+    if (debug)
+    {
+        Pout<< "label pointBoundaryMesh::findPatchID(const word&) const"
+            << " Patch named " << patchName << " not found.  "
+            << "Available patch names: " << names() << endl;
+    }
+
+    // Not found, return -1
+    return -1;
 }
 
 
@@ -243,4 +702,80 @@ void Foam::pointBoundaryMesh::updateMesh()
 }
 
 
+void Foam::pointBoundaryMesh::reorder
+(
+    const labelUList& oldToNew,
+    const bool validBoundary
+)
+{
+    // Change order of patches
+    pointPatchList::reorder(oldToNew);
+
+    // Adapt indices
+    pointPatchList& patches = *this;
+
+    forAll(patches, patchi)
+    {
+        patches[patchi].index() = patchi;
+    }
+
+    // Clear group-to-patch addressing. Note: could re-calculate
+    groupIDsPtr_.reset(nullptr);
+
+    if (validBoundary)
+    {
+        updateMesh();
+    }
+
+    if (debug)
+    {
+        pointPatchList& Patches = *this;
+
+        Pout<< "pointBoundaryMesh::reorder"
+            << "(const labelUList&, const bool): "
+            << "reordered pointBoundaryMesh:" << endl;
+        Pout<< incrIndent;
+        for (const auto& pp : Patches)
+        {
+            Pout<< indent
+                << "index:" << pp.index() << " patch:" << pp.name()
+                << " type:" << pp.type() << endl;
+        }
+        Pout<< decrIndent;
+    }
+}
+
+
+bool Foam::pointBoundaryMesh::writeData(Ostream& os) const
+{
+    const pointPatchList& patches = *this;
+
+    os  << patches.size() << nl << token::BEGIN_LIST << incrIndent << nl;
+
+    forAll(patches, patchi)
+    {
+        os  << indent << patches[patchi].name() << nl
+            << indent << token::BEGIN_BLOCK << nl
+            << incrIndent << patches[patchi] << decrIndent
+            << indent << token::END_BLOCK << endl;
+    }
+
+    os  << decrIndent << token::END_LIST;
+
+    // Check state of IOstream
+    os.check("pointBoundaryMesh::writeData(Ostream& os) const");
+
+    return os.good();
+}
+
+
+//bool Foam::pointBoundaryMesh::writeObject
+//(
+//    IOstreamOption
+//) const
+//{
+//    return regIOobject::writeObject(fmt, ver, IOstream::UNCOMPRESSED);
+//}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/pointMesh/pointBoundaryMesh/pointBoundaryMesh.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2018-2023 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,6 +39,7 @@ SourceFiles
 #define Foam_pointBoundaryMesh_H
 
 #include "pointPatch.H"
+#include "regIOobject.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -56,19 +57,29 @@ class wordRes;
 
 class pointBoundaryMesh
 :
-    public pointPatchList
+    public pointPatchList,
+    public regIOobject
 {
     // Private Data
 
         //- Reference to mesh
         const pointMesh& mesh_;
 
+        //- Demand-driven: list of patch ids per group
+        mutable autoPtr<HashTable<labelList>> groupIDsPtr_;
+
 
     // Private Member Functions
 
         //- Calculate geometry for the patches (transformation tensors etc.)
         void calcGeometry();
 
+        //- Some patches have inGroup entries
+        bool hasGroupIDs() const;
+
+        //- Calculate group name to patch ids lookup
+        void calcGroupIDs() const;
+
         //- Assign facePointPatches corresponding to the given polyBoundaryMesh
         void addPatches(const polyBoundaryMesh& pbm);
 
@@ -85,11 +96,27 @@ public:
     friend class pointMesh;
 
 
+    //- Runtime type information
+    TypeName("pointBoundaryMesh");
+
+
     // Constructors
 
         //- Construct from polyBoundaryMesh
         pointBoundaryMesh(const pointMesh&, const polyBoundaryMesh&);
 
+        //- Construct from IOobject and polyBoundaryMesh
+        pointBoundaryMesh
+        (
+            const IOobject& io,
+            const pointMesh&,
+            const polyBoundaryMesh&
+        );
+
+
+    //- Destructor
+    virtual ~pointBoundaryMesh() = default;
+
 
     // Member Functions
 
@@ -99,6 +126,21 @@ public:
             return mesh_;
         }
 
+        //- The number of patches before the first processor patch.
+        label nNonProcessor() const;
+
+        //- The number of processorPointPatch patches
+        label nProcessorPatches() const;
+
+        //- Return a list of patch names
+        wordList names() const;
+
+        //- Return a list of patch types
+        wordList types() const;
+
+        //- Return a list of physical types
+        wordList physicalTypes() const;
+
         //- Return (sorted) patch indices for all matches.
         //  A no-op (returns empty list) for an empty matcher
         labelList indices(const wordRe& matcher, const bool useGroups) const;
@@ -121,14 +163,38 @@ public:
 
         //- Find patch index given a name
         //  A no-op (returns -1) for an empty patchName
-        label findPatchID(const word& patchName) const;
+        label findPatchID
+        (
+            const word& patchName,
+            const bool allowNotFound = true
+        ) const;
 
-        //- Correct polyBoundaryMesh after moving points
+        //- The patch indices per patch group
+        const HashTable<labelList>& groupPatchIDs() const;
+
+        //- Correct pointBoundaryMesh after moving points
         void movePoints(const pointField&);
 
-        //- Correct polyBoundaryMesh after topology update
+        //- Correct pointBoundaryMesh after topology update
         void updateMesh();
 
+        //- Reorders patches. Ordering does not have to be done in
+        //  ascending or descending order. Reordering has to be unique.
+        //  (is shuffle) If validBoundary calls updateMesh()
+        //  after reordering to recalculate data (so call needs to be parallel
+        //  sync in that case)
+        void reorder(const labelUList& oldToNew, const bool validBoundary);
+
+        //- writeData member function required by regIOobject
+        virtual bool writeData(Ostream&) const;
+
+//        //- Write using given format, version and form uncompression
+//        virtual bool writeObject
+//        (
+//            IOstreamOption streamOpt,
+//            const bool writeOnProc
+//        ) const;
+
 
     // Housekeeping
 

src/OpenFOAM/meshes/pointMesh/pointMesh.C

@@ -37,9 +37,12 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(pointMesh, 0);
+    defineTypeNameAndDebug(pointMesh, 0);
 }
 
+Foam::word Foam::pointMesh::meshSubDir = "pointMesh";
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::pointMesh::mapFields(const mapPolyMesh& mpm)
@@ -88,8 +91,66 @@ Foam::pointMesh::pointMesh(const polyMesh& pMesh)
 }
 
 
+Foam::pointMesh::pointMesh(const polyMesh& pMesh, const IOobject& io)
+:
+    MeshObject<polyMesh, Foam::UpdateableMeshObject, pointMesh>(pMesh),
+    GeoMesh<polyMesh>(pMesh),
+    boundary_(io, *this, pMesh.boundaryMesh())
+{
+    if (debug)
+    {
+        Pout<< "pointMesh::pointMesh(const polyMesh&): "
+            << "Constructing from IO " << io.objectRelPath()
+            << endl;
+    }
+
+    // Calculate the geometry for the patches (transformation tensors etc.)
+    boundary_.calcGeometry();
+}
+
+
+Foam::pointMesh::pointMesh
+(
+    const polyMesh& pMesh,
+    const IOobjectOption::readOption rOpt
+)
+:
+    pointMesh
+    (
+        pMesh,
+        IOobject
+        (
+            pMesh.name(),                // polyMesh region
+            pMesh.facesInstance(),       // polyMesh topology instance
+            pMesh.time(),
+            rOpt,
+            Foam::IOobject::NO_WRITE
+        )
+    )
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+void Foam::pointMesh::setInstance
+(
+    const fileName& inst,
+    const IOobjectOption::writeOption wOpt
+)
+{
+    if (debug)
+    {
+        Pout<< "pointMesh::setInstance(): "
+            << "Setting instance to " << inst << endl;
+    }
+    this->writeOpt(wOpt);
+    this->instance() = inst;
+
+    boundary_.writeOpt(wOpt);
+    boundary_.instance() = inst;
+}
+
+
 bool Foam::pointMesh::movePoints()
 {
     if (debug)
@@ -118,4 +179,19 @@ void Foam::pointMesh::updateMesh(const mapPolyMesh& mpm)
 }
 
 
+bool Foam::pointMesh::writeObject
+(
+    IOstreamOption streamOpt,
+    const bool writeOnProc
+) const
+{
+    if (debug)
+    {
+        Pout<< "pointMesh::writeObject(IOstreamOption, const bool): "
+            << "Writing to " << boundary_.objectRelPath() << endl;
+    }
+    return boundary_.writeObject(streamOpt, writeOnProc);
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/pointMesh/pointMesh.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
-    Copyright (C) 2021-2023 OpenCFD Ltd.
+    Copyright (C) 2021-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -86,12 +86,26 @@ public:
     // Declare name of the class and its debug switch
     ClassName("pointMesh");
 
+    //- Return the mesh sub-directory name (usually "pointMesh")
+    static word meshSubDir;
+
 
     // Constructors
 
         //- Construct from polyMesh
         explicit pointMesh(const polyMesh& pMesh);
 
+        //- Construct from polyMesh and IOobject (used when reading boundary)
+        explicit pointMesh(const polyMesh& pMesh, const IOobject& io);
+
+        //- Construct from polyMesh and readOpt. Takes instance, time etc
+        //- from polyMesh. Used when reading boundary.
+        explicit pointMesh
+        (
+            const polyMesh& pMesh,
+            const IOobjectOption::readOption rOpt
+        );
+
 
     //- Destructor
     ~pointMesh() = default;
@@ -141,6 +155,13 @@ public:
             return GeoMesh<polyMesh>::mesh_.time();
         }
 
+        //- Set the instance for mesh files
+        void setInstance
+        (
+            const fileName& instance,
+            const IOobjectOption::writeOption wOpt = IOobject::AUTO_WRITE
+        );
+
 
     // Volume Mesh
 
@@ -171,6 +192,16 @@ public:
         {
             return &pm == this;
         }
+
+
+    // Write
+
+        //- Write
+        virtual bool writeObject
+        (
+            IOstreamOption streamOpt,
+            const bool writeOnProc = true
+        ) const;
 };
 
 

src/OpenFOAM/meshes/pointMesh/pointPatches/basic/generic/genericPointPatch.H

@@ -72,6 +72,39 @@ public:
         :
             facePointPatch(patch, bm)
         {}
+
+        //- Construct given the original patch and a map
+        genericPointPatch
+        (
+            const genericPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        )
+        :
+            facePointPatch(patch, bm, index, mapAddressing, reversePointMap)
+        {}
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<genericPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
 };
 
 

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.C

@@ -91,6 +91,20 @@ Foam::cyclicPointPatch::cyclicPointPatch
 {}
 
 
+Foam::cyclicPointPatch::cyclicPointPatch
+(
+    const cyclicPointPatch& patch,
+    const pointBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const labelUList& reversePointMap
+)
+:
+    coupledFacePointPatch(patch, bm, index, mapAddressing, reversePointMap),
+    cyclicPolyPatch_(refCast<const cyclicPolyPatch>(patch.patch()))
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::cyclicPointPatch::~cyclicPointPatch()

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclic/cyclicPointPatch.H

@@ -105,6 +105,36 @@ public:
             const pointBoundaryMesh& bm
         );
 
+        //- Construct given the original patch and a map
+        cyclicPointPatch
+        (
+            const cyclicPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        );
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<cyclicPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     //- Destructor
     virtual ~cyclicPointPatch();

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/cyclicSlip/cyclicSlipPointPatch.H

@@ -73,6 +73,39 @@ public:
             cyclicPointPatch(patch, bm)
         {}
 
+        //- Construct given the original patch and a map
+        cyclicSlipPointPatch
+        (
+            const cyclicSlipPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        )
+        :
+            cyclicPointPatch(patch, bm, index, mapAddressing, reversePointMap)
+        {}
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<cyclicSlipPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     //- Destructor
     virtual ~cyclicSlipPointPatch() = default;

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/empty/emptyPointPatch.H

@@ -72,6 +72,39 @@ public:
             facePointPatch(patch, bm)
         {}
 
+        //- Construct given the original patch and a map
+        emptyPointPatch
+        (
+            const emptyPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        )
+        :
+            facePointPatch(patch, bm, index, mapAddressing, reversePointMap)
+        {}
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<emptyPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     // Member Functions
 

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/nonuniformTransformCyclic/nonuniformTransformCyclicPointPatch.H

@@ -73,6 +73,40 @@ public:
             cyclicPointPatch(patch, bm)
         {}
 
+        //- Construct given the original patch and a map
+        nonuniformTransformCyclicPointPatch
+        (
+            const nonuniformTransformCyclicPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        )
+        :
+            cyclicPointPatch(patch, bm, index, mapAddressing, reversePointMap)
+        {}
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::
+                NewFrom<nonuniformTransformCyclicPointPatch>
+                (
+                    *this,
+                    bm,
+                    index,
+                    mapAddressing,
+                    reversePointMap
+                );
+        }
+
 
     // Destructor
 

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.C

@@ -117,4 +117,20 @@ Foam::processorPointPatch::processorPointPatch
 {}
 
 
+Foam::processorPointPatch::processorPointPatch
+(
+    const processorPointPatch& patch,
+    const pointBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const labelUList& reversePointMap
+)
+:
+    coupledFacePointPatch(patch, bm, index, mapAddressing, reversePointMap),
+    procPolyPatch_(refCast<const processorPolyPatch>(patch.patch()))
+{
+    //? map reverseMeshPoints_ or leave demand-driven
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processor/processorPointPatch.H

@@ -111,6 +111,36 @@ public:
             const pointBoundaryMesh& bm
         );
 
+        //- Construct given the original patch and a map
+        processorPointPatch
+        (
+            const processorPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        );
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<processorPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     //- Destructor
     virtual ~processorPointPatch() = default;

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processorCyclic/processorCyclicPointPatch.C

@@ -59,6 +59,20 @@ processorCyclicPointPatch::processorCyclicPointPatch
 {}
 
 
+Foam::processorCyclicPointPatch::processorCyclicPointPatch
+(
+    const processorCyclicPointPatch& patch,
+    const pointBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const labelUList& reversePointMap
+)
+:
+    processorPointPatch(patch, bm, index, mapAddressing, reversePointMap),
+    procCycPolyPatch_(refCast<const processorCyclicPolyPatch>(patch.patch()))
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 processorCyclicPointPatch::~processorCyclicPointPatch()

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/processorCyclic/processorCyclicPointPatch.H

@@ -89,6 +89,36 @@ public:
             const pointBoundaryMesh& bm
         );
 
+        //- Construct given the original patch and a map
+        processorCyclicPointPatch
+        (
+            const processorCyclicPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        );
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<processorCyclicPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     //- Destructor
     virtual ~processorCyclicPointPatch();

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/symmetry/symmetryPointPatch.H

@@ -75,6 +75,38 @@ public:
             facePointPatch(patch, bm)
         {}
 
+        //- Construct given the original patch and a map
+        symmetryPointPatch
+        (
+            const symmetryPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ):
+            facePointPatch(patch, bm, index, mapAddressing, reversePointMap)
+        {}
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<symmetryPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     // Member Functions
 

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/symmetryPlane/symmetryPlanePointPatch.C

@@ -58,6 +58,23 @@ Foam::symmetryPlanePointPatch::symmetryPlanePointPatch
 {}
 
 
+Foam::symmetryPlanePointPatch::symmetryPlanePointPatch
+(
+    const symmetryPlanePointPatch& patch,
+    const pointBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const labelUList& reversePointMap
+)
+:
+    facePointPatch(patch, bm, index, mapAddressing, reversePointMap),
+    symmetryPlanePolyPatch_
+    (
+        refCast<const symmetryPlanePolyPatch>(patch.patch())
+    )
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::symmetryPlanePointPatch::applyConstraint

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/symmetryPlane/symmetryPlanePointPatch.H

@@ -74,6 +74,36 @@ public:
             const pointBoundaryMesh& bm
         );
 
+        //- Construct given the original patch and a map
+        symmetryPlanePointPatch
+        (
+            const symmetryPlanePointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        );
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<symmetryPlanePointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     // Member Functions
 

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.C

@@ -58,6 +58,20 @@ Foam::wedgePointPatch::wedgePointPatch
 {}
 
 
+Foam::wedgePointPatch::wedgePointPatch
+(
+    const wedgePointPatch& patch,
+    const pointBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const labelUList& reversePointMap
+)
+:
+    facePointPatch(patch, bm, index, mapAddressing, reversePointMap),
+    wedgePolyPatch_(refCast<const wedgePolyPatch>(patch.patch()))
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::wedgePointPatch::applyConstraint

src/OpenFOAM/meshes/pointMesh/pointPatches/constraint/wedge/wedgePointPatch.H

@@ -74,6 +74,36 @@ public:
             const pointBoundaryMesh& bm
         );
 
+        //- Construct given the original patch and a map
+        wedgePointPatch
+        (
+            const wedgePointPatch& pp,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        );
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<wedgePointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     // Member Functions
 

src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.C

@@ -50,4 +50,19 @@ Foam::coupledFacePointPatch::coupledFacePointPatch
 {}
 
 
+Foam::coupledFacePointPatch::coupledFacePointPatch
+(
+    const coupledFacePointPatch& patch,
+    const pointBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const labelUList& reversePointMap
+)
+:
+    facePointPatch(patch, bm, index, mapAddressing, reversePointMap),
+    coupledPointPatch(bm),
+    coupledPolyPatch_(refCast<const coupledPolyPatch>(patch.patch()))
+{}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/pointMesh/pointPatches/derived/coupled/coupledFacePointPatch.H

@@ -91,6 +91,16 @@ public:
             const pointBoundaryMesh& bm
         );
 
+        //- Construct given the original patch and a map
+        coupledFacePointPatch
+        (
+            const coupledFacePointPatch& pp,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        );
+
 
     //- Destructor
     virtual ~coupledFacePointPatch() = default;

src/OpenFOAM/meshes/pointMesh/pointPatches/derived/wall/wallPointPatch.H

@@ -71,6 +71,39 @@ public:
         :
             facePointPatch(patch, bm)
         {}
+
+        //- Construct given the original patch and a map
+        wallPointPatch
+        (
+            const wallPointPatch& patch,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        )
+        :
+            facePointPatch(patch, bm, index, mapAddressing, reversePointMap)
+        {}
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<wallPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
 };
 
 

src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.C

@@ -84,9 +84,30 @@ Foam::facePointPatch::facePointPatch
     const pointBoundaryMesh& bm
 )
 :
-    pointPatch(bm),
+    pointPatch(p.name(), p.index(), bm, p.physicalType(), p.inGroups()),
     polyPatch_(p)
 {}
 
 
+Foam::facePointPatch::facePointPatch
+(
+    const facePointPatch& pp,
+    const pointBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const labelUList& reversePointMap
+)
+:
+    pointPatch
+    (
+        pp.name(),
+        index,
+        bm,
+        pp.physicalType(),
+        pp.inGroups()
+    ),
+    polyPatch_(pp.patch())
+{}
+
+
 // ************************************************************************* //

src/OpenFOAM/meshes/pointMesh/pointPatches/facePointPatch/facePointPatch.H

@@ -128,6 +128,36 @@ public:
             const pointBoundaryMesh& pm
         );
 
+        //- Construct given the original patch and a map
+        facePointPatch
+        (
+            const facePointPatch& pp,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        );
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<facePointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
 
     // Selectors
 

src/OpenFOAM/meshes/pointMesh/pointPatches/meshPointPatch/meshPointPatch.C

@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "meshPointPatch.H"
+#include "addToRunTimeSelectionTable.H"
+#include "pointMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(meshPointPatch, 0);
+
+//- Needs run-time selection table on pointPatch, not facePointPatch
+addToRunTimeSelectionTable
+(
+    pointPatch,
+    meshPointPatch,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::meshPointPatch::meshPointPatch
+(
+    const word& name,
+    const labelUList& meshPoints,
+    const vectorField& pointNormals,
+    const label index,
+    const pointBoundaryMesh& bm,
+    const word& patchType
+)
+:
+    pointPatch(name, index, bm, word::null, wordList()),
+    meshPoints_(meshPoints),
+    pointNormals_(pointNormals)
+{
+    if (meshPoints_.size() != pointNormals_.size())
+    {
+        FatalErrorInFunction << "patch " << name
+            << " size of meshPoints " << meshPoints_.size()
+            << " differs from size of pointNormals " << pointNormals_.size()
+            << exit(FatalError);
+    }
+}
+
+
+Foam::meshPointPatch::meshPointPatch
+(
+    const word& name,
+    const dictionary& dict,
+    const label index,
+    const pointBoundaryMesh& bm,
+    const word& patchType
+)
+:
+    pointPatch(name, dict, index, bm),
+    meshPoints_(dict.get<labelList>("meshPoints")),
+    pointNormals_("normals", dict, meshPoints_.size())
+{}
+
+
+Foam::meshPointPatch::meshPointPatch
+(
+    const meshPointPatch& pp,
+    const pointBoundaryMesh& bm,
+    const label index,
+    const labelUList& mapAddressing,
+    const labelUList& reversePointMap
+)
+:
+    meshPointPatch
+    (
+        pp.name(),
+        labelList(reversePointMap, labelList(pp.meshPoints(), mapAddressing)),
+        vectorField(pp.pointNormals(), mapAddressing),
+        index,
+        bm,
+        pp.type()
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::meshPointPatch::movePoints(PstreamBuffers&, const pointField& p)
+{
+    localPointsPtr_.reset(nullptr);
+
+    // Recalculate the point normals? Something like
+    //if (owner())
+    //{
+    //    const SubList<face> subFaces
+    //    (
+    //        mesh.faces(),
+    //        mesh.nBoundaryFaces(),
+    //        mesh.nInternalFaces()
+    //    );
+    //    const primitivePatch pp(subFaces, mesh.points());
+    //
+    //
+    //    for (const label pointi : meshPoints())
+    //    {
+    //        const auto fnd(pp.meshPointMap().find(pointi));
+    //        if (fnd)
+    //        {
+    //            const label patchPointi = fnd();
+    //            // Determine point patch equiv
+    //
+    //        const auto& point
+    //
+    //
+
+}
+    
+
+void Foam::meshPointPatch::updateMesh(PstreamBuffers&)
+{
+    localPointsPtr_.reset(nullptr);
+    // Do what to pointNormals? Don't know what the new mesh points are
+}
+
+
+const Foam::pointField& Foam::meshPointPatch::localPoints() const
+{
+    if (!localPointsPtr_)
+    {
+        localPointsPtr_.reset
+        (
+            new pointField
+            (
+                boundaryMesh().mesh().mesh().points(),
+                meshPoints()
+            )
+        );
+    }
+    return localPointsPtr_();
+}
+
+
+//const Foam::vectorField& Foam::meshPointPatch::pointNormals() const
+//{
+//    if (!pointNormalsPtr_)
+//    {
+//        pointNormalsPtr_.reset(new vectorField(size()));
+//        vectorField& pointNormals = pointNormalsPtr_();
+//        forAll(constraints_, i)
+//        {
+//            pointNormals[i] = constraints_[i].second();
+//        }
+//    }
+//    return pointNormalsPtr_();
+//}
+
+
+//void Foam::meshPointPatch::applyConstraint
+//(
+//    const label pointi,
+//    pointConstraint& pc
+//) const
+//{
+//    pc.combine(constraints_[pointi]);
+//}
+//
+//
+//void Foam::meshPointPatch::setConstraints
+//(
+//    const List<pointConstraint>& pc
+//)
+//{
+//    constraints_ = pc;
+//    localPointsPtr_.reset(nullptr);
+//    pointNormalsPtr_.reset(nullptr);
+//}
+
+
+void Foam::meshPointPatch::write(Ostream& os) const
+{
+    pointPatch::write(os);
+    meshPoints().writeEntry("meshPoints", os);
+    pointNormals().writeEntry("normals", os);
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/meshes/pointMesh/pointPatches/meshPointPatch/meshPointPatch.H

@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2016 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::meshPointPatch
+
+Description
+    pointPatch with explicitly provided points instead of using the points
+    of a polyPatch.
+
+SourceFiles
+    meshPointPatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef meshPointPatch_H
+#define meshPointPatch_H
+
+#include "pointPatch.H"
+#include "polyPatch.H"
+#include "autoPtr.H"
+#include "patchIdentifier.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                  Class meshPointPatch Declaration
+\*---------------------------------------------------------------------------*/
+
+class meshPointPatch
+:
+    public pointPatch
+{
+private:
+
+    // Private Member Functions
+
+        //- No copy construct
+        meshPointPatch(const meshPointPatch&) = delete;
+
+        //- No copy assignment
+        void operator=(const meshPointPatch&) = delete;
+
+
+protected:
+
+    // Protected Member Data
+
+        //- Explicit mesh points
+        const labelList meshPoints_;
+
+        const vectorField pointNormals_;
+
+        //- Demand-driven local points
+        mutable autoPtr<pointField> localPointsPtr_;
+
+
+    // Protected Member Functions
+
+        //- Correct patches after moving points
+        virtual void movePoints(PstreamBuffers&, const pointField&);
+
+        //- Update of the patch topology
+        virtual void updateMesh(PstreamBuffers&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("meshPoint");
+
+
+    // Constructors
+
+        //- Construct from components
+        meshPointPatch
+        (
+            const word& name,
+            const labelUList& meshPoints,
+            const vectorField& pointNormals,
+            const label index,
+            const pointBoundaryMesh& bm,
+            const word& patchType
+        );
+
+        //- Construct from dictionary
+        meshPointPatch
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const pointBoundaryMesh& bm,
+            const word& patchType
+        );
+
+        //- Construct given the original patch and a map
+        meshPointPatch
+        (
+            const meshPointPatch& pp,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        );
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const
+        {
+            return autoPtr<pointPatch>::NewFrom<meshPointPatch>
+            (
+                *this,
+                bm,
+                index,
+                mapAddressing,
+                reversePointMap
+            );
+        }
+
+
+    //- Destructor
+    virtual ~meshPointPatch() = default;
+
+
+    // Member Functions
+
+        //- Return size
+        virtual label size() const
+        {
+            return meshPoints().size();
+        }
+
+        //- Return mesh points
+        virtual const labelList& meshPoints() const
+        {
+            return meshPoints_;
+        }
+
+        //- Return pointField of points in patch
+        virtual const pointField& localPoints() const;
+
+        //- Return point unit normals
+        virtual const vectorField& pointNormals() const
+        {
+            return pointNormals_;
+        }
+
+        //- Write the pointPatch data as a dictionary
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/meshes/pointMesh/pointPatches/pointPatch/pointPatch.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,7 +32,33 @@ License
 
 namespace Foam
 {
-defineTypeNameAndDebug(pointPatch, 0);
+    defineTypeNameAndDebug(pointPatch, 0);
+
+//    int pointPatch::disallowGenericPointPatch
+//    (
+//        debug::debugSwitch("disallowGenericPointPatch", 0)
+//    );
+
+    defineRunTimeSelectionTable(pointPatch, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::pointPatch::write(Ostream& os) const
+{
+    os.writeKeyword("type") << type() << token::END_STATEMENT << nl;
+    patchIdentifier::write(os);
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const pointPatch& p)
+{
+    p.write(os);
+    os.check("Ostream& operator<<(Ostream& os, const pointPatch& p");
+    return os;
 }
 
 

src/OpenFOAM/meshes/pointMesh/pointPatches/pointPatch/pointPatch.H

@@ -38,6 +38,7 @@ SourceFiles
 #ifndef Foam_pointPatch_H
 #define Foam_pointPatch_H
 
+#include "patchIdentifier.H"
 #include "labelList.H"
 #include "vectorField.H"
 #include "pointField.H"
@@ -54,6 +55,8 @@ class pointConstraint;
 class pointPatch;
 class PstreamBuffers;
 
+Ostream& operator<<(Ostream&, const pointPatch&);
+
 //- Store lists of pointPatch as a PtrList
 typedef PtrList<pointPatch> pointPatchList;
 
@@ -62,6 +65,8 @@ typedef PtrList<pointPatch> pointPatchList;
 \*---------------------------------------------------------------------------*/
 
 class pointPatch
+:
+    public patchIdentifier
 {
     // Private Data
 
@@ -114,14 +119,90 @@ public:
     TypeName("basePatch");
 
 
+    // Declare run-time constructor selection tables
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            pointPatch,
+            dictionary,
+            (
+                const word& name,
+                const dictionary& dict,
+                const label index,
+                const pointBoundaryMesh& bm,
+                const word& patchType
+            ),
+            (name, dict, index, bm, patchType)
+        );
+
+
     // Constructor
 
-        //- Construct from boundary mesh
-        explicit pointPatch(const pointBoundaryMesh& bm)
+        //- Construct from components
+        explicit pointPatch
+        (
+            const word& name,
+            const label index,
+            const pointBoundaryMesh& bm,
+            const word& physicalType,
+            const wordList& inGroups
+        )
         :
+            patchIdentifier(name, index, physicalType, inGroups),
             boundaryMesh_(bm)
         {}
 
+        //- Construct from dictionary
+        explicit pointPatch
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const pointBoundaryMesh& bm
+        )
+        :
+            patchIdentifier(name, dict, index),
+            boundaryMesh_(bm)
+        {}
+
+        //- Construct given the original patch and a map
+        explicit pointPatch
+        (
+            const pointPatch& pp,
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        )
+        :
+            patchIdentifier(pp.name(), index, pp.physicalType(), pp.inGroups()),
+            boundaryMesh_(bm)
+        {}
+
+        //- Construct and return a subset clone,
+        //- resetting the point list and boundary mesh
+        virtual autoPtr<pointPatch> clone
+        (
+            const pointBoundaryMesh& bm,
+            const label index,
+            const labelUList& mapAddressing,
+            const labelUList& reversePointMap
+        ) const = 0;
+
+
+    // Selectors
+
+        //- Return a pointer to a new patch created on freestore. Returns
+        //- null if not found
+        static autoPtr<pointPatch> New
+        (
+            const word& name,
+            const dictionary& dict,
+            const label index,
+            const pointBoundaryMesh&
+        );
+
 
     //- Destructor
     virtual ~pointPatch() = default;
@@ -129,15 +210,9 @@ public:
 
     // Member Functions
 
-        //- Return name
-        virtual const word& name() const = 0;
-
         //- Return size
         virtual label size() const = 0;
 
-        //- Return the index of this patch in the pointBoundaryMesh
-        virtual label index() const = 0;
-
         //- Return boundaryMesh reference
         const pointBoundaryMesh& boundaryMesh() const
         {
@@ -172,6 +247,14 @@ public:
             pointConstraint&
         ) const
         {}
+
+        //- Write the pointPatch data as a dictionary
+        virtual void write(Ostream&) const;
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<<(Ostream&, const pointPatch&);
 };
 
 

src/OpenFOAM/meshes/pointMesh/pointPatches/pointPatch/pointPatchNew.C

@@ -0,0 +1,56 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2016 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::pointPatch> Foam::pointPatch::New
+(
+    const word& name,
+    const dictionary& dict,
+    const label index,
+    const pointBoundaryMesh& bm
+)
+{
+    // Similar to polyPatchNew but no support for generic since we want it
+    // to fall through to the construct-from-polyPatch
+    DebugInFunction << "Constructing pointPatch" << endl;
+
+    const word patchType(dict.lookup("type"));
+    //dict.readIfPresent("geometricType", patchType);
+
+    auto* ctorPtr = dictionaryConstructorTable(patchType);
+
+    if (!ctorPtr)
+    {
+        return nullptr;
+    }
+
+    return autoPtr<pointPatch>(ctorPtr(name, dict, index, bm, patchType));
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/meshes/polyMesh/polyBoundaryMesh/polyBoundaryMesh.C

@@ -1404,6 +1404,9 @@ void Foam::polyBoundaryMesh::reorder
         patches[patchi].index() = patchi;
     }
 
+    // Clear group-to-patch addressing. Note: could re-calculate
+    groupIDsPtr_.reset(nullptr);
+
     if (validBoundary)
     {
         updateMesh();

src/OpenFOAM/primitives/strings/lists/SubStrings.H

@@ -109,19 +109,14 @@ public:
             this->push_back(range);
         }
 
-        //- Const reference to the first element,
-        //- for consistency with other OpenFOAM containers
-        auto first() const -> decltype(this->front())
-        {
-            return this->front();
-        }
-
-        //- Const reference to the last element,
-        //- for consistency with other OpenFOAM containers
-        auto last() const -> decltype(this->back())
-        {
-            return this->back();
-        }
+    // FUTURE?
+    //  #if __cplusplus >= 201703L
+    //  std::string_view view(size_t pos) const
+    //  {}
+    //  #else
+    //  stdFoam::span<const char> view(size_t pos) const
+    //  {}
+    //  #endif
 };
 
 

src/dynamicMesh/motionSolvers/displacement/displacement/displacementMotionSolver.C

@@ -57,7 +57,7 @@ Foam::displacementMotionSolver::displacementMotionSolver
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        pointMesh::New(mesh)
+        pointMesh::New(mesh, Foam::IOobject::READ_IF_PRESENT)
     )
 {}
 

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/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2022,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,10 +32,16 @@ License
 #include "cyclicPolyPatch.H"
 #include "emptyPolyPatch.H"
 #include "processorPolyPatch.H"
+#include "meshPointPatch.H"
+#include "processorPointPatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-Foam::word Foam::fvMeshSubset::exposedPatchName("oldInternalFaces");
+namespace Foam
+{
+    word fvMeshSubset::exposedPatchName("oldInternalFaces");
+    defineTypeNameAndDebug(fvMeshSubset, 0);
+}
 
 
 // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
@@ -522,6 +528,16 @@ void Foam::fvMeshSubset::reset
 
 void Foam::fvMeshSubset::reset(const Foam::zero)
 {
+    // Was old pointMesh present?
+    const auto* basePointMeshPtr =
+            baseMesh_.thisDb().cfindObject<pointMesh>(pointMesh::typeName);
+    if (basePointMeshPtr)
+    {
+        DebugPout<< "fvMeshSubset::reset(const Foam::zero) :"
+                << " Detected pointMesh" << endl;
+    }
+
+
     clear();
 
     // Create zero-sized subMesh
@@ -573,6 +589,46 @@ void Foam::fvMeshSubset::reset(const Foam::zero)
     }
 
 
+    // Clone old additional point patches
+    if (basePointMeshPtr)
+    {
+        DebugPout<< "Subsetting pointMesh" << endl;
+        const auto& basePointMesh = *basePointMeshPtr;
+        const auto& oldPointBoundary = basePointMesh.boundary();
+
+        // 1. Generate pointBoundaryMesh from polyBoundaryMesh (so ignoring
+        //    any additional patches
+        const auto& newSubPointMesh = pointMesh::New(newSubMesh);
+
+        auto& newBoundary =
+            const_cast<pointBoundaryMesh&>(newSubPointMesh.boundary());
+
+        // Start off from (poly)patch map
+        pointPatchMap_ = patchMap_;
+
+        // 2. Explicitly add subsetted meshPointPatches
+        for (const auto& oldPointPatch : oldPointBoundary)
+        {
+            const auto* mppPtr = isA<meshPointPatch>(oldPointPatch);
+            if (mppPtr && (newBoundary.findPatchID(mppPtr->name()) == -1))
+            {
+                newBoundary.push_back
+                (
+                    mppPtr->clone
+                    (
+                        newBoundary,
+                        newBoundary.size(),
+                        labelList::null(), // map
+                        labelList::null()  // map
+                    )
+                );
+            }
+        }
+
+        // Extend patchMap with -1
+        pointPatchMap_.setSize(newBoundary.size(), -1);
+    }
+
     // Add the zones
     subsetZones();
 }
@@ -585,6 +641,16 @@ void Foam::fvMeshSubset::reset
     const bool syncPar
 )
 {
+    // Was old pointMesh present?
+    const auto* basePointMeshPtr =
+            baseMesh_.thisDb().cfindObject<pointMesh>(pointMesh::typeName);
+    if (basePointMeshPtr)
+    {
+        DebugPout<< "fvMeshSubset::reset(const bitSet&) :"
+                << " Detected pointMesh" << endl;
+    }
+
+
     // Clear all old maps and pointers
     clear();
 
@@ -1125,6 +1191,8 @@ void Foam::fvMeshSubset::reset
 
     // Inserted patch
 
+    label newInternalPatchID = -1;
+
     if (wantedPatchID == -1)
     {
         label oldInternalSize = boundaryPatchSizes[oldInternalPatchID];
@@ -1158,6 +1226,7 @@ void Foam::fvMeshSubset::reset
             // the internal faces
             patchStart += boundaryPatchSizes[oldInternalPatchID];
             patchMap_[nNewPatches] = -1;
+            newInternalPatchID = nNewPatches;
             ++nNewPatches;
         }
     }
@@ -1232,6 +1301,98 @@ void Foam::fvMeshSubset::reset
 
     // Subset and add any zones
     subsetZones();
+
+
+    if (basePointMeshPtr)
+    {
+        DebugPout<< "Subsetting pointMesh" << endl;
+        const auto& basePointMesh = *basePointMeshPtr;
+        const auto& oldPointBoundary = basePointMesh.boundary();
+
+        // 1. Generate pointBoundaryMesh from polyBoundaryMesh (so ignoring
+        //    any additional patches
+        const auto& newSubPointMesh = pointMesh::New(subMeshPtr_());
+
+        pointPatchMap_ = patchMap_;
+
+        auto& newBoundary =
+            const_cast<pointBoundaryMesh&>(newSubPointMesh.boundary());
+
+
+        // 2. Explicitly add subsetted meshPointPatches
+        labelList oldToNewPoints(baseMesh_.nPoints(), -1);
+        forAll(pointMap_, i)
+        {
+            oldToNewPoints[pointMap_[i]] = i;
+        }
+
+
+        // Add meshPointPatches
+        pointPatchMap_.setSize(newBoundary.size(), -1);
+
+        for (const auto& oldPointPatch : oldPointBoundary)
+        {
+            const auto* mppPtr = isA<meshPointPatch>(oldPointPatch);
+            if (mppPtr && (newBoundary.findPatchID(mppPtr->name()) == -1))
+            {
+                const auto& mp = mppPtr->meshPoints();
+                DynamicList<label> subPointMap(mp.size());
+                forAll(mp, i)
+                {
+                    const label newPointi = oldToNewPoints[mp[i]];
+                    if (newPointi != -1)
+                    {
+                        subPointMap.append(i);
+                    }
+                }
+
+                pointPatchMap_.push_back(mppPtr->index());
+
+                newBoundary.push_back
+                (
+                    mppPtr->clone
+                    (
+                        newBoundary,
+                        newBoundary.size(),
+                        subPointMap,        // map
+                        oldToNewPoints
+                    )
+                );
+            }
+        }
+
+
+        // 3. rotate into place:
+        //      - global patches (including meshPointPatches)
+        //      - optional 'internalFaces' patch
+        //      - processor patches
+        labelList oldToNew(newBoundary.size());
+        label newPatchi = 0;
+        forAll(newBoundary, patchi)
+        {
+            if
+            (
+                patchi != newInternalPatchID
+            && !isA<processorPointPatch>(newBoundary[patchi])
+            )
+            {
+                oldToNew[patchi] = newPatchi++;
+            }
+        }
+        if (newInternalPatchID != -1)
+        {
+            oldToNew[newInternalPatchID] = newPatchi++;
+        }
+        forAll(newBoundary, patchi)
+        {
+            if (isA<processorPointPatch>(newBoundary[patchi]))
+            {
+                oldToNew[patchi] = newPatchi++;
+            }
+        }
+        newBoundary.reorder(oldToNew, true);
+        inplaceReorder(oldToNew, pointPatchMap_);
+    }
 }
 
 

src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.H

@@ -102,6 +102,9 @@ class fvMeshSubset
         //- Patch mapping array
         labelList patchMap_;
 
+        //- PointPatch mapping array
+        labelList pointPatchMap_;
+
 
     // Private Member Functions
 
@@ -135,6 +138,10 @@ protected:
 
 public:
 
+    // Declare name of the class and its debug switch
+    ClassName("fvMeshSubset");
+
+
     // Static Data Members
 
         //- Name for exposed internal faces (default: oldInternalFaces)
@@ -225,6 +232,10 @@ public:
         //- Return patch map
         inline const labelList& patchMap() const;
 
+        //- Return point-patch map. Usually identical to patchMap except if
+        //- additional patches are added to the pointMesh.
+        inline const labelList& pointPatchMap() const;
+
 
     // Edit
 

src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetI.H

@@ -104,4 +104,19 @@ inline const Foam::labelList& Foam::fvMeshSubset::patchMap() const
 }
 
 
+inline const Foam::labelList& Foam::fvMeshSubset::pointPatchMap() const
+{
+    checkHasSubMesh();
+
+    if (pointPatchMap_.empty())
+    {
+        return patchMap_;
+    }
+    else
+    {
+        return pointPatchMap_;
+    }
+}
+
+
 // ************************************************************************* //

src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubsetTemplates.C

@@ -554,7 +554,7 @@ Foam::fvMeshSubset::interpolate
     (
         sf,
         pointMesh::New(subMesh()),     // subsetted point mesh
-        patchMap(),
+        pointPatchMap(),
         pointMap()
     );
 }

src/functionObjects/field/wallHeatFlux/wallHeatFlux.C

@@ -108,6 +108,10 @@ Foam::functionObjects::wallHeatFlux::wallHeatFlux
     writeFile(obr_, name, typeName, dict),
     qrName_("qr")
 {
+    read(dict);
+
+    writeFileHeader(file());
+
     volScalarField* wallHeatFluxPtr
     (
         new volScalarField
@@ -127,10 +131,6 @@ Foam::functionObjects::wallHeatFlux::wallHeatFlux
     );
 
     mesh_.objectRegistry::store(wallHeatFluxPtr);
-
-    read(dict);
-
-    writeFileHeader(file());
 }
 
 

src/fvMotionSolver/Make/files

@@ -38,6 +38,8 @@ $(derivedPoint)/angularOscillatingDisplacement/angularOscillatingDisplacementPoi
 $(derivedPoint)/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
 $(derivedPoint)/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
 $(derivedPoint)/waveDisplacement/waveDisplacementPointPatchVectorField.C
+$(derivedPoint)/edgeSlipDisplacement/edgeSlipDisplacementPointPatchVectorField.C
+$(derivedPoint)/pointAttractionDisplacement/pointAttractionDisplacementPointPatchVectorField.C
 
 $(derivedPoint)/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchFields.C
 $(derivedPoint)/uniformInterpolatedDisplacement/uniformInterpolatedDisplacementPointPatchVectorField.C

src/fvMotionSolver/fvMotionSolvers/componentDisplacement/componentLaplacian/displacementComponentLaplacianFvMotionSolver.C

@@ -145,6 +145,18 @@ Foam::displacementComponentLaplacianFvMotionSolver::curPoints() const
         pointDisplacement_
     );
 
+    // Evaluate the bcs so they are consistent with the internal field
+    // Might fight the multi-patch behaviour inside volPointInterpolate
+    if
+    (
+        pointDisplacement_.boundaryField().size()
+     != cellDisplacement_.boundaryField().size()
+    )
+    {
+        pointDisplacement_.correctBoundaryConditions();
+    }
+
+
     if (pointLocation_)
     {
         if (debug)

src/fvMotionSolver/fvMotionSolvers/displacement/SBRStress/displacementSBRStressFvMotionSolver.C

@@ -160,6 +160,17 @@ Foam::displacementSBRStressFvMotionSolver::curPoints() const
         pointDisplacement_
     );
 
+    // Evaluate the bcs so they are consistent with the internal field
+    // Might fight the multi-patch behaviour inside volPointInterpolate
+    if
+    (
+        pointDisplacement_.boundaryField().size()
+     != cellDisplacement_.boundaryField().size()
+    )
+    {
+        pointDisplacement_.correctBoundaryConditions();
+    }
+
     tmp<pointField> tcurPoints
     (
         points0() + pointDisplacement().primitiveField()
@@ -186,6 +197,10 @@ void Foam::displacementSBRStressFvMotionSolver::solve()
        *diffusivityPtr_->operator()()
     );
 
+    // Make sure the cellMotion bcs are consistent with the pointDisplacement.
+    // This makes sure the grad below uses more up-to-date values.
+    cellDisplacement_.boundaryFieldRef().updateCoeffs();
+
     volTensorField gradCd("gradCd", fvc::grad(cellDisplacement_));
 
     fv::options& fvOptions(fv::options::New(fvMesh_));

src/fvMotionSolver/fvMotionSolvers/displacement/laplacian/displacementLaplacianFvMotionSolver.C

@@ -267,6 +267,17 @@ Foam::displacementLaplacianFvMotionSolver::curPoints() const
         pointDisplacement_
     );
 
+    // Evaluate the bcs so they are consistent with the internal field
+    // Might fight the multi-patch behaviour inside volPointInterpolate
+    if
+    (
+        pointDisplacement_.boundaryField().size()
+     != cellDisplacement_.boundaryField().size()
+    )
+    {
+        pointDisplacement_.correctBoundaryConditions();
+    }
+
     if (pointLocation_)
     {
         if (debug)
@@ -333,6 +344,10 @@ void Foam::displacementLaplacianFvMotionSolver::solve()
     diffusivity().correct();
     pointDisplacement_.boundaryFieldRef().updateCoeffs();
 
+    // Make sure the cellMotion bcs are consistent with the pointDisplacement.
+    // This makes sure the grad below uses more up-to-date values.
+    cellDisplacement_.boundaryFieldRef().updateCoeffs();
+
     fv::options& fvOptions(fv::options::New(fvMesh_));
 
     // We explicitly do NOT want to interpolate the motion inbetween

src/fvMotionSolver/fvMotionSolvers/displacement/solidBodyDisplacementLaplacian/solidBodyDisplacementLaplacianFvMotionSolver.C

@@ -271,6 +271,17 @@ Foam::solidBodyDisplacementLaplacianFvMotionSolver::curPoints() const
         pointDisplacement_
     );
 
+    // Evaluate the bcs so they are consistent with the internal field
+    // Might fight the multi-patch behaviour inside volPointInterpolate
+    if
+    (
+        pointDisplacement_.boundaryField().size()
+     != cellDisplacement_.boundaryField().size()
+    )
+    {
+        pointDisplacement_.correctBoundaryConditions();
+    }
+
     tmp<pointField> tnewPoints
     (
         transformPoints(SBMFPtr_().transformation(), points0())

src/fvMotionSolver/fvMotionSolvers/fvMotionSolver/fvMotionSolverTemplates.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
+    Copyright (C) 2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,6 +29,7 @@ License
 #include "fvMotionSolver.H"
 #include "fixedValuePointPatchFields.H"
 #include "cellMotionFvPatchFields.H"
+#include "facePointPatch.H"
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
@@ -38,23 +40,33 @@ Foam::wordList Foam::fvMotionSolver::cellMotionBoundaryTypes
     Boundary& pmUbf
 ) const
 {
-    wordList cmUbf = pmUbf.types();
+    wordList cmUbf(fvMesh_.boundary().size());
 
-    // Remove global patches from the end of the list
-    cmUbf.setSize(fvMesh_.boundary().size());
-
-    forAll(cmUbf, patchi)
+    forAll(pmUbf, patchi)
     {
-        if (isA<fixedValuePointPatchField<Type>>(pmUbf[patchi]))
+        const auto& pfld = pmUbf[patchi];
+        const auto* fppPtr = isA<facePointPatch>(pfld.patch());
+        if (fppPtr)
         {
-            cmUbf[patchi] = cellMotionFvPatchField<Type>::typeName;
-        }
+            const auto& fpp = *fppPtr;
+            const label polyPatchi = fpp.patch().index();
 
-        if (debug)
-        {
-            Pout<< "Patch:" << fvMesh_.boundary()[patchi].patch().name()
-                << " pointType:" << pmUbf.types()[patchi]
-                << " cellType:" << cmUbf[patchi] << endl;
+            if (isA<fixedValuePointPatchField<Type>>(pfld))
+            {
+                cmUbf[polyPatchi] = cellMotionFvPatchField<Type>::typeName;
+            }
+            else
+            {
+                // Take over pointPatch type
+                cmUbf[polyPatchi] = pfld.type();
+            }
+
+            if (debug)
+            {
+                Pout<< "Patch:" << fvMesh_.boundary()[patchi].patch().name()
+                    << " pointType:" << pfld.type()
+                    << " cellType:" << cmUbf[patchi] << endl;
+            }
         }
     }
 

src/fvMotionSolver/fvPatchFields/derived/cellMotion/cellMotionFvPatchField.H

@@ -37,7 +37,6 @@ SourceFiles
 #ifndef cellMotionFvPatchField_H
 #define cellMotionFvPatchField_H
 
-#include "Random.H"
 #include "fixedValueFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/fvMotionSolver/pointPatchFields/derived/edgeSlipDisplacement/edgeSlipDisplacementPointPatchVectorField.C

@@ -0,0 +1,391 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "edgeSlipDisplacementPointPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "transformField.H"
+#include "displacementMotionSolver.H"
+#include "featureEdgeMesh.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::edgeSlipDisplacementPointPatchVectorField::read
+(
+    const objectRegistry& obr,
+    const dictionary& dict
+)
+{
+    const Time& tm = obr.time();
+
+    const fileName featFileName(dict.get<fileName>("file", keyType::LITERAL));
+
+    if (tm.foundObject<edgeMesh>(featFileName))
+    {
+        return;
+    }
+
+    IOobject extFeatObj
+    (
+        featFileName,                       // name
+        tm.constant(),                      // instance
+        "extendedFeatureEdgeMesh",          // local
+        obr,                                // registry
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        IOobject::REGISTER
+    );
+
+    //const fileName fName(typeFilePath<extendedFeatureEdgeMesh>(extFeatObj));
+    const fileName fName(extFeatObj.typeFilePath<extendedFeatureEdgeMesh>());
+
+    if (!fName.empty() && extendedEdgeMesh::canRead(fName))
+    {
+        Info<< "Reading edgeMesh from " << extFeatObj.objectRelPath() << endl;
+        auto* eMeshPtr = new extendedFeatureEdgeMesh(extFeatObj);
+        eMeshPtr->store();
+    }
+    else
+    {
+        // Try reading as edgeMesh
+
+        IOobject featObj
+        (
+            featFileName,                       // name
+            tm.constant(),                      // instance
+            "triSurface",                       // local
+            obr,                                // registry
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            IOobject::REGISTER
+        );
+
+        Info<< "Reading edgeMesh from " << featObj.objectRelPath() << endl;
+        const fileName fName(featObj.typeFilePath<featureEdgeMesh>());
+
+        if (fName.empty())
+        {
+            FatalIOErrorInFunction(dict)
+                << "Could not open " << featObj.objectPath()
+                << exit(FatalIOError);
+        }
+
+        // Read as edgeMesh
+        auto* eMeshPtr = new featureEdgeMesh(featObj);
+        eMeshPtr->store();
+    }
+}
+
+
+void Foam::edgeSlipDisplacementPointPatchVectorField::calcProjection
+(
+    vectorField& displacement
+) const
+{
+    const polyMesh& mesh = patch().boundaryMesh().mesh()();
+    const pointField& localPoints = patch().localPoints();
+    const labelList& meshPoints = patch().meshPoints();
+
+    //const scalar deltaT = mesh.time().deltaTValue();
+
+    // Construct large enough vector in direction of projectDir so
+    // we're guaranteed to hit something.
+
+    //- Per point projection vector:
+    const scalar projectLen = mesh.bounds().mag();
+
+
+
+    // Get fixed points (bit of a hack)
+    const pointZone* zonePtr = nullptr;
+
+    if (frozenPointsZone_.size() > 0)
+    {
+        const pointZoneMesh& pZones = mesh.pointZones();
+
+        zonePtr = &pZones[frozenPointsZone_];
+
+        Info<< "edgeSlipDisplacementPointPatchVectorField : Fixing all "
+            << zonePtr->size() << " points in pointZone " << zonePtr->name()
+            << endl;
+    }
+
+    // Get the starting locations from the motionSolver
+    const pointField& points0 = mesh.lookupObject<displacementMotionSolver>
+    (
+        "dynamicMeshDict"
+    ).points0();
+
+
+    pointField start(meshPoints.size());
+    forAll(start, i)
+    {
+        start[i] = points0[meshPoints[i]] + displacement[i];
+    }
+
+    const auto& tree = edgeTree();
+
+    label nNotProjected = 0;
+    forAll(meshPoints, i)
+    {
+        if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
+        {
+            // Fixed point. Reset to point0 location.
+            displacement[i] = points0[meshPoints[i]] - localPoints[i];
+        }
+        else
+        {
+            pointIndexHit nearest = tree.findNearest(start[i], sqr(projectLen));
+            if (nearest.hit())
+            {
+                displacement[i] = nearest.point() - points0[meshPoints[i]];
+            }
+            else
+            {
+                nNotProjected++;
+
+                if (debug)
+                {
+                    Pout<< "    point:" << meshPoints[i]
+                        << " coord:" << localPoints[i]
+                        << "  did not find any surface within " << projectLen
+                        << endl;
+                }
+            }
+        }
+    }
+
+    reduce(nNotProjected, sumOp<label>());
+
+    if (nNotProjected > 0)
+    {
+        Info<< "edgeSlipDisplacement :"
+            << " on patch " << patch().name()
+            << " did not project " << nNotProjected
+            << " out of " << returnReduce(localPoints.size(), sumOp<label>())
+            << " points." << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::edgeSlipDisplacementPointPatchVectorField::
+edgeSlipDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    pointPatchVectorField(p, iF),
+    velocity_(Zero)
+{}
+
+
+Foam::edgeSlipDisplacementPointPatchVectorField::
+edgeSlipDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const dictionary& dict
+)
+:
+    pointPatchVectorField(p, iF, dict),
+    velocity_(dict.get<vector>("velocity")),
+    featFileName_(dict.get<fileName>("file", keyType::LITERAL)),
+    frozenPointsZone_(dict.getOrDefault("frozenPointsZone", word::null))
+{
+    read(this->patch().boundaryMesh().mesh().time(), dict);
+}
+
+
+Foam::edgeSlipDisplacementPointPatchVectorField::
+edgeSlipDisplacementPointPatchVectorField
+(
+    const edgeSlipDisplacementPointPatchVectorField& ppf,
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const pointPatchFieldMapper&
+)
+:
+    pointPatchVectorField(p, iF),
+    velocity_(ppf.velocity_),
+    featFileName_(ppf.featFileName_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+Foam::edgeSlipDisplacementPointPatchVectorField::
+edgeSlipDisplacementPointPatchVectorField
+(
+    const edgeSlipDisplacementPointPatchVectorField& ppf
+)
+:
+    pointPatchVectorField(ppf),
+    velocity_(ppf.velocity_),
+    featFileName_(ppf.featFileName_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+Foam::edgeSlipDisplacementPointPatchVectorField::
+edgeSlipDisplacementPointPatchVectorField
+(
+    const edgeSlipDisplacementPointPatchVectorField& ppf,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    pointPatchVectorField(ppf, iF),
+    velocity_(ppf.velocity_),
+    featFileName_(ppf.featFileName_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::indexedOctree<Foam::treeDataEdge>&
+Foam::edgeSlipDisplacementPointPatchVectorField::edgeTree() const
+{
+    if (!edgeTreePtr_)
+    {
+        const Time& tm = this->patch().boundaryMesh().mesh().time();
+        const auto& eMesh = tm.lookupObject<edgeMesh>(featFileName_);
+
+        const pointField& points = eMesh.points();
+        const edgeList& edges = eMesh.edges();
+
+        // Calculate bb of all points
+        treeBoundBox bb(points);
+
+        // Random number generator. Bit dodgy since not exactly random ;-)
+        Random rndGen(65431);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        bb.inflate(rndGen, 1e-4, ROOTVSMALL);
+
+        edgeTreePtr_.reset
+        (
+            new indexedOctree<treeDataEdge>
+            (
+                treeDataEdge(edges, points),  // All edges
+
+                bb,     // overall search domain
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+    }
+
+    return edgeTreePtr_();
+}
+
+
+void Foam::edgeSlipDisplacementPointPatchVectorField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const vectorField currentDisplacement(this->patchInternalField());
+
+    // Calculate displacement to project points onto surface
+    vectorField displacement(currentDisplacement);
+    calcProjection(displacement);
+
+    // offset wrt current displacement
+    displacement -= currentDisplacement;
+
+    // Clip offset to maximum displacement possible: velocity*timestep
+
+    const Time& tm = this->patch().boundaryMesh().mesh().time();
+    const scalar deltaT = tm.deltaTValue();
+    const vector clipVelocity = velocity_*deltaT;
+
+    forAll(displacement, i)
+    {
+        vector& d = displacement[i];
+
+        for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
+        {
+            if (d[cmpt] < 0)
+            {
+                d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
+            }
+            else
+            {
+                d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
+            }
+        }
+    }
+
+
+    // Get internal field to insert values into
+    Field<vector>& iF = const_cast<Field<vector>&>(this->primitiveField());
+
+    displacement += currentDisplacement;
+
+    setInInternalField(iF, displacement);
+
+    pointPatchVectorField::updateCoeffs();
+}
+
+
+void Foam::edgeSlipDisplacementPointPatchVectorField::write
+(
+    Ostream& os
+) const
+{
+    pointPatchField<vector>::write(os);
+    os.writeEntry("file", featFileName_);
+    os.writeEntryIfDifferent<word>
+    (
+        "frozenPointsZone",
+        word::null,
+        frozenPointsZone_
+    );
+    os.writeEntry("velocity", velocity_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+makePointPatchTypeField
+(
+    pointPatchVectorField,
+    edgeSlipDisplacementPointPatchVectorField
+);
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/fvMotionSolver/pointPatchFields/derived/edgeSlipDisplacement/edgeSlipDisplacementPointPatchVectorField.H

@@ -0,0 +1,198 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::edgeSlipDisplacementPointPatchVectorField
+
+Description
+    Displacement follows an edgeMesh. Use in a displacementMotionSolver
+    as a bc on the pointDisplacement field.
+
+    Needs:
+    - file : name of edgeMesh. Either:
+        - .extendedFeatureEdgeMesh (read from  constant/extendedFeatureEdgeMesh)
+        - or .eMesh in or a known format (read from constant/triSurface)
+    - frozenPointsZone : empty or name of pointZone containing points
+                         that do not move
+    - velocity : maximum displacement velocity
+
+
+SourceFiles
+    edgeSlipDisplacementPointPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef edgeSlipDisplacementPointPatchVectorField_H
+#define edgeSlipDisplacementPointPatchVectorField_H
+
+#include "pointPatchFields.H"
+#include "extendedFeatureEdgeMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+          Class edgeSlipDisplacementPointPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class edgeSlipDisplacementPointPatchVectorField
+:
+    public pointPatchVectorField
+{
+private:
+
+    // Private data
+
+        //- Maximum velocity
+        const vector velocity_;
+
+        //- file
+        const fileName featFileName_;
+
+        //- pointZone with frozen points
+        const word frozenPointsZone_;
+
+        //- Edge searching
+        mutable autoPtr<indexedOctree<treeDataEdge>> edgeTreePtr_;
+
+
+    // Private Member Functions
+
+        //- Calculate displacement (w.r.t. points0()) to project onto surface
+        void calcProjection(vectorField& displacement) const;
+
+        //- No copy assignment
+        void operator=
+        (
+            const edgeSlipDisplacementPointPatchVectorField&
+        ) = delete;
+
+
+protected:
+
+        const indexedOctree<treeDataEdge>& edgeTree() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("edgeSlipDisplacement");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        edgeSlipDisplacementPointPatchVectorField
+        (
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        edgeSlipDisplacementPointPatchVectorField
+        (
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given patch field onto a new patch
+        edgeSlipDisplacementPointPatchVectorField
+        (
+            const edgeSlipDisplacementPointPatchVectorField&,
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&,
+            const pointPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        edgeSlipDisplacementPointPatchVectorField
+        (
+            const edgeSlipDisplacementPointPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<pointPatchVectorField> clone() const
+        {
+            return autoPtr<pointPatchVectorField>
+            (
+                new edgeSlipDisplacementPointPatchVectorField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        edgeSlipDisplacementPointPatchVectorField
+        (
+            const edgeSlipDisplacementPointPatchVectorField&,
+            const DimensionedField<vector, pointMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual autoPtr<pointPatchVectorField> clone
+        (
+            const DimensionedField<vector, pointMesh>& iF
+        ) const
+        {
+            return autoPtr<pointPatchVectorField>
+            (
+                new edgeSlipDisplacementPointPatchVectorField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+    // Member Functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Read (& store) geometry. Exposed so point attraction can reuse it.
+        static void read
+        (
+            const objectRegistry& obr,
+            const dictionary& dict
+        );
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvMotionSolver/pointPatchFields/derived/pointAttractionDisplacement/pointAttractionDisplacementPointPatchVectorField.C

@@ -0,0 +1,329 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointAttractionDisplacementPointPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "transformField.H"
+#include "displacementMotionSolver.H"
+#include "featureEdgeMesh.H"
+#include "edgeSlipDisplacementPointPatchVectorField.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::pointAttractionDisplacementPointPatchVectorField::calcProjection
+(
+    vectorField& displacement
+) const
+{
+    const polyMesh& mesh = patch().boundaryMesh().mesh()();
+    const pointField& localPoints = patch().localPoints();
+    const labelList& meshPoints = patch().meshPoints();
+
+    //const scalar deltaT = mesh.time().deltaTValue();
+
+    // Construct large enough vector in direction of projectDir so
+    // we're guaranteed to hit something.
+
+    //- Per point projection vector:
+    const scalar projectLen = mesh.bounds().mag();
+
+
+
+    // Get fixed points (bit of a hack)
+    const pointZone* zonePtr = nullptr;
+
+    if (frozenPointsZone_.size() > 0)
+    {
+        const pointZoneMesh& pZones = mesh.pointZones();
+
+        zonePtr = &pZones[frozenPointsZone_];
+
+        Pout<< "pointAttractionDisplacementPointPatchVectorField : Fixing all "
+            << zonePtr->size() << " points in pointZone " << zonePtr->name()
+            << endl;
+    }
+
+    // Get the starting locations from the motionSolver
+    const pointField& points0 = mesh.lookupObject<displacementMotionSolver>
+    (
+        "dynamicMeshDict"
+    ).points0();
+
+
+    pointField start(meshPoints.size());
+    forAll(start, i)
+    {
+        start[i] = points0[meshPoints[i]] + displacement[i];
+    }
+
+    const auto& tree = pointTree();
+
+    label nNotProjected = 0;
+    forAll(meshPoints, i)
+    {
+        if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
+        {
+            // Fixed point. Reset to point0 location.
+            displacement[i] = points0[meshPoints[i]] - localPoints[i];
+        }
+        else
+        {
+            pointIndexHit nearest = tree.findNearest(start[i], sqr(projectLen));
+            if (nearest.hit())
+            {
+                displacement[i] = nearest.point() - points0[meshPoints[i]];
+            }
+            else
+            {
+                nNotProjected++;
+
+                if (debug)
+                {
+                    Pout<< "    point:" << meshPoints[i]
+                        << " coord:" << localPoints[i]
+                        << "  did not find any surface within " << projectLen
+                        << endl;
+                }
+            }
+        }
+    }
+
+    reduce(nNotProjected, sumOp<label>());
+
+    if (nNotProjected > 0)
+    {
+        Info<< "pointAttractionDisplacement :"
+            << " on patch " << patch().name()
+            << " did not project " << nNotProjected
+            << " out of " << returnReduce(localPoints.size(), sumOp<label>())
+            << " points." << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pointAttractionDisplacementPointPatchVectorField::
+pointAttractionDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    pointPatchVectorField(p, iF),
+    velocity_(Zero)
+{}
+
+
+Foam::pointAttractionDisplacementPointPatchVectorField::
+pointAttractionDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const dictionary& dict
+)
+:
+    pointPatchVectorField(p, iF, dict),
+    velocity_(dict.get<vector>("velocity")),
+    featFileName_(dict.get<fileName>("file", keyType::LITERAL)),
+    frozenPointsZone_(dict.getOrDefault("frozenPointsZone", word::null))
+{
+    // Read&store edge mesh on registry
+    edgeSlipDisplacementPointPatchVectorField::read
+    (
+        this->patch().boundaryMesh().mesh().time(),
+        dict
+    );
+}
+
+
+Foam::pointAttractionDisplacementPointPatchVectorField::
+pointAttractionDisplacementPointPatchVectorField
+(
+    const pointAttractionDisplacementPointPatchVectorField& ppf,
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const pointPatchFieldMapper&
+)
+:
+    pointPatchVectorField(p, iF),
+    velocity_(ppf.velocity_),
+    featFileName_(ppf.featFileName_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+Foam::pointAttractionDisplacementPointPatchVectorField::
+pointAttractionDisplacementPointPatchVectorField
+(
+    const pointAttractionDisplacementPointPatchVectorField& ppf
+)
+:
+    pointPatchVectorField(ppf),
+    velocity_(ppf.velocity_),
+    featFileName_(ppf.featFileName_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+Foam::pointAttractionDisplacementPointPatchVectorField::
+pointAttractionDisplacementPointPatchVectorField
+(
+    const pointAttractionDisplacementPointPatchVectorField& ppf,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    pointPatchVectorField(ppf, iF),
+    velocity_(ppf.velocity_),
+    featFileName_(ppf.featFileName_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::indexedOctree<Foam::treeDataPoint>&
+Foam::pointAttractionDisplacementPointPatchVectorField::pointTree() const
+{
+    if (!pointTreePtr_)
+    {
+        const Time& tm = this->patch().boundaryMesh().mesh().time();
+        const auto& eMesh = tm.lookupObject<edgeMesh>(featFileName_);
+
+        const pointField& points = eMesh.points();
+
+        // Calculate bb of all points
+        treeBoundBox bb(points);
+
+        // Random number generator. Bit dodgy since not exactly random ;-)
+        Random rndGen(65431);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        bb.inflate(rndGen, 1e-4, ROOTVSMALL);
+
+        pointTreePtr_.reset
+        (
+            new indexedOctree<treeDataPoint>
+            (
+                treeDataPoint(points),  // All edges
+
+                bb,     // overall search domain
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+        );
+    }
+
+    return pointTreePtr_();
+}
+
+
+void Foam::pointAttractionDisplacementPointPatchVectorField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const vectorField currentDisplacement(this->patchInternalField());
+
+    // Calculate displacement to project points onto surface
+    vectorField displacement(currentDisplacement);
+    calcProjection(displacement);
+
+    // offset wrt current displacement
+    displacement -= currentDisplacement;
+
+    // Clip offset to maximum displacement possible: velocity*timestep
+
+    const Time& tm = this->patch().boundaryMesh().mesh().time();
+    const scalar deltaT = tm.deltaTValue();
+    const vector clipVelocity = velocity_*deltaT;
+
+    forAll(displacement, i)
+    {
+        vector& d = displacement[i];
+
+        for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
+        {
+            if (d[cmpt] < 0)
+            {
+                d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
+            }
+            else
+            {
+                d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
+            }
+        }
+    }
+
+
+    // Get internal field to insert values into
+    Field<vector>& iF = const_cast<Field<vector>&>(this->primitiveField());
+
+    displacement += currentDisplacement;
+
+    setInInternalField(iF, displacement);
+
+    pointPatchVectorField::updateCoeffs();
+}
+
+
+void Foam::pointAttractionDisplacementPointPatchVectorField::write
+(
+    Ostream& os
+) const
+{
+    pointPatchField<vector>::write(os);
+    os.writeEntry("file", featFileName_);
+    os.writeEntryIfDifferent<word>
+    (
+        "frozenPointsZone",
+        word::null,
+        frozenPointsZone_
+    );
+    os.writeEntry("velocity", velocity_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+makePointPatchTypeField
+(
+    pointPatchVectorField,
+    pointAttractionDisplacementPointPatchVectorField
+);
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/fvMotionSolver/pointPatchFields/derived/pointAttractionDisplacement/pointAttractionDisplacementPointPatchVectorField.H

@@ -0,0 +1,191 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2024 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::pointAttractionDisplacementPointPatchVectorField
+
+Description
+    Displacement by attraction to nearest point. Use in a
+    displacementMotionSolver as a bc on the pointDisplacement field.
+
+    Needs:
+    - file : name of edgeMesh. Either:
+        - .extendedFeatureEdgeMesh (read from  constant/extendedFeatureEdgeMesh)
+        - or .eMesh in or a known format (read from constant/triSurface)
+    - frozenPointsZone : empty or name of pointZone containing points
+                         that do not move
+    - velocity : maximum displacement velocity
+
+
+SourceFiles
+    pointAttractionDisplacementPointPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointAttractionDisplacementPointPatchVectorField_H
+#define pointAttractionDisplacementPointPatchVectorField_H
+
+#include "pointPatchFields.H"
+#include "extendedFeatureEdgeMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+      Class pointAttractionDisplacementPointPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointAttractionDisplacementPointPatchVectorField
+:
+    public pointPatchVectorField
+{
+private:
+
+    // Private data
+
+        //- Maximum velocity
+        const vector velocity_;
+
+        //- file
+        const fileName featFileName_;
+
+        //- pointZone with frozen points
+        const word frozenPointsZone_;
+
+        //- Point searching
+        mutable autoPtr<indexedOctree<treeDataPoint>> pointTreePtr_;
+
+
+    // Private Member Functions
+
+        //- Calculate displacement (w.r.t. points0()) to project onto surface
+        void calcProjection(vectorField& displacement) const;
+
+        //- No copy assignment
+        void operator=
+        (
+            const pointAttractionDisplacementPointPatchVectorField&
+        ) = delete;
+
+
+protected:
+
+        const indexedOctree<treeDataPoint>& pointTree() const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("pointAttraction");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        pointAttractionDisplacementPointPatchVectorField
+        (
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        pointAttractionDisplacementPointPatchVectorField
+        (
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given patch field onto a new patch
+        pointAttractionDisplacementPointPatchVectorField
+        (
+            const pointAttractionDisplacementPointPatchVectorField&,
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&,
+            const pointPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        pointAttractionDisplacementPointPatchVectorField
+        (
+            const pointAttractionDisplacementPointPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<pointPatchVectorField> clone() const
+        {
+            return autoPtr<pointPatchVectorField>
+            (
+                new pointAttractionDisplacementPointPatchVectorField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        pointAttractionDisplacementPointPatchVectorField
+        (
+            const pointAttractionDisplacementPointPatchVectorField&,
+            const DimensionedField<vector, pointMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual autoPtr<pointPatchVectorField> clone
+        (
+            const DimensionedField<vector, pointMesh>& iF
+        ) const
+        {
+            return autoPtr<pointPatchVectorField>
+            (
+                new pointAttractionDisplacementPointPatchVectorField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+    // Member Functions
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C

@@ -318,7 +318,12 @@ surfaceDisplacementPointPatchVectorField
     velocity_(dict.get<vector>("velocity")),
     surfacesDict_(dict.subDict("geometry")),
     projectMode_(projectModeNames_.get("projectMode", dict)),
-    projectDir_(dict.get<vector>("projectDirection")),
+    projectDir_
+    (
+        (projectMode_ == FIXEDNORMAL)
+      ? dict.get<vector>("projectDirection")
+      : Zero
+    ),
     wedgePlane_(dict.getOrDefault("wedgePlane", -1)),
     frozenPointsZone_(dict.getOrDefault("frozenPointsZone", word::null))
 {
@@ -424,12 +429,13 @@ void Foam::surfaceDisplacementPointPatchVectorField::updateCoeffs()
 
     const polyMesh& mesh = patch().boundaryMesh().mesh()();
 
-    vectorField currentDisplacement(this->patchInternalField());
+    const vectorField currentDisplacement(this->patchInternalField());
 
     // Calculate intersections with surface w.r.t points0.
     vectorField displacement(currentDisplacement);
     calcProjection(displacement);
 
+
     // offset wrt current displacement
     vectorField offset(displacement-currentDisplacement);
 
@@ -442,21 +448,15 @@ void Foam::surfaceDisplacementPointPatchVectorField::updateCoeffs()
     {
         vector& d = offset[i];
 
-        for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
+        const scalar magD(mag(d));
+        if (magD > ROOTVSMALL)
         {
-            if (d[cmpt] < 0)
-            {
-                d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
-            }
-            else
-            {
-                d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
-            }
+            d /= magD;
+            d *= min(magD, mag(clipVelocity));
         }
     }
 
     this->operator==(currentDisplacement+offset);
-
     fixedValuePointPatchVectorField::updateCoeffs();
 }
 
@@ -467,9 +467,13 @@ void Foam::surfaceDisplacementPointPatchVectorField::write(Ostream& os) const
     os.writeEntry("velocity", velocity_);
     os.writeEntry("geometry", surfacesDict_);
     os.writeEntry("projectMode", projectModeNames_[projectMode_]);
-    os.writeEntry("projectDirection", projectDir_);
-    os.writeEntry("wedgePlane", wedgePlane_);
-
+    os.writeEntryIfDifferent<vector>
+    (
+        "projectDirection",
+        Zero,
+        projectDir_
+    );
+    os.writeEntryIfDifferent<label>("wedgePlane", -1, wedgePlane_);
     os.writeEntryIfDifferent<word>
     (
         "frozenPointsZone",

src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C

@@ -312,7 +312,12 @@ surfaceSlipDisplacementPointPatchVectorField
     pointPatchVectorField(p, iF, dict),
     surfacesDict_(dict.subDict("geometry")),
     projectMode_(projectModeNames_.get("projectMode", dict)),
-    projectDir_(dict.get<vector>("projectDirection")),
+    projectDir_
+    (
+        (projectMode_ == FIXEDNORMAL)
+      ? dict.get<vector>("projectDirection")
+      : Zero
+    ),
     wedgePlane_(dict.getOrDefault("wedgePlane", -1)),
     frozenPointsZone_(dict.getOrDefault("frozenPointsZone", word::null))
 {}
@@ -425,9 +430,13 @@ void Foam::surfaceSlipDisplacementPointPatchVectorField::write
     pointPatchField<vector>::write(os);
     os.writeEntry("geometry", surfacesDict_);
     os.writeEntry("projectMode", projectModeNames_[projectMode_]);
-    os.writeEntry("projectDirection", projectDir_);
-    os.writeEntry("wedgePlane", wedgePlane_);
-
+    os.writeEntryIfDifferent<vector>
+    (
+        "projectDirection",
+        Zero,
+        projectDir_
+    );
+    os.writeEntryIfDifferent<label>("wedgePlane", -1, wedgePlane_);
     os.writeEntryIfDifferent<word>
     (
         "frozenPointsZone",

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.C

@@ -90,6 +90,8 @@ void Foam::fv::actuationDiskSource::writeFileHeader(Ostream& os)
         writeFile::writeCommented(os, "P");
     }
 
+    writeFile::writeCommented(os, "diskDir");
+
     os  << endl;
 }
 
@@ -208,14 +210,7 @@ Foam::fv::actuationDiskSource::actuationDiskSource
             scalarMinMax::ge(VSMALL)
         )
     ),
-    diskDir_
-    (
-        coeffs_.getCheck<vector>
-        (
-            "diskDir",
-            [&](const vector& vec){ return mag(vec) > VSMALL; }
-        ).normalise()
-    ),
+    diskDir_(Function1<vector>::New("diskDir", coeffs_, &mesh)),
     UvsCpPtr_(Function1<scalar>::New("Cp", coeffs_, &mesh)),
     UvsCtPtr_(Function1<scalar>::New("Ct", coeffs_, &mesh)),
     monitorCells_()
@@ -237,6 +232,25 @@ Foam::fv::actuationDiskSource::actuationDiskSource
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::vector Foam::fv::actuationDiskSource::diskDir() const
+{
+    const scalar t = mesh_.time().timeOutputValue();
+    const vector dir(diskDir_->value(t));
+    const scalar magDir = Foam::mag(dir);
+
+    if (magDir < SMALL)
+    {
+        FatalErrorInFunction
+            << "Actuator disk surface-normal vector is zero: " << dir
+            << " at time=" << t
+            << exit(FatalError);
+    }
+
+    // normalised:
+    return dir/magDir;
+}
+
+
 void Foam::fv::actuationDiskSource::addSup
 (
     fvMatrix<vector>& eqn,
@@ -295,19 +309,10 @@ bool Foam::fv::actuationDiskSource::read(const dictionary& dict)
                 << exit(FatalIOError);
         }
 
-        dict.readIfPresent("diskDir", diskDir_);
-        diskDir_.normalise();
-        if (mag(diskDir_) < VSMALL)
-        {
-            FatalIOErrorInFunction(dict)
-                << "Actuator disk surface-normal vector is zero: "
-                << "diskDir = " << diskDir_
-                << exit(FatalIOError);
-        }
+        // TBD: runTime re-reading of "diskDir" ?
 
         return true;
     }
-
     return false;
 }
 

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSource.H

@@ -170,7 +170,7 @@ Usage
       type       | Type name: actuationDiskSource      | word | yes      | -
       diskArea   | Actuator disk planar surface area   | scalar | yes    | -
       diskDir    | Surface-normal vector of the actuator disk <!--
-               --> pointing downstream                 | vector    | yes | -
+               --> pointing downstream                 | Function1 | yes | -
       Cp         | Power coefficient                   | Function1 | yes | -
       Ct         | Thrust coefficient                  | Function1 | yes | -
       monitorMethod | Type of incoming velocity monitoring method   <!--
@@ -207,8 +207,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef actuationDiskSource_H
-#define actuationDiskSource_H
+#ifndef Foam_actuationDiskSource_H
+#define Foam_actuationDiskSource_H
 
 #include "cellSetOption.H"
 #include "writeFile.H"
@@ -276,7 +276,8 @@ protected:
         scalar diskArea_;
 
         //- Surface-normal vector of the actuator disk pointing downstream
-        vector diskDir_;
+        //- as a function of time
+        autoPtr<Function1<vector>> diskDir_;
 
         //- Velocity vs power coefficients
         autoPtr<Function1<scalar>> UvsCpPtr_;
@@ -364,14 +365,11 @@ public:
 
         // Access
 
-            //- Normal disk direction
-            const vector& diskDir() const
-            {
-                return diskDir_;
-            }
+            //- Normal disk direction, evaluated at timeOutputValue
+            vector diskDir() const;
 
             //- Disk area
-            scalar diskArea() const
+            scalar diskArea() const noexcept
             {
                 return diskArea_;
             }

src/fvOptions/sources/derived/actuationDiskSource/actuationDiskSourceTemplates.C

@@ -110,12 +110,13 @@ void Foam::fv::actuationDiskSource::calcFroudeMethod
     }
 
     // (BJSB:Eq. 3.9)
+    const vector diskDir = this->diskDir();
     const scalar a = 1.0 - Cp/Ct;
-    const scalar T = 2.0*rhoRef*diskArea_*magSqr(Uref & diskDir_)*a*(1 - a);
+    const scalar T = 2.0*rhoRef*diskArea_*magSqr(Uref & diskDir)*a*(1 - a);
 
     for (const label celli : cells_)
     {
-        Usource[celli] += ((cellsV[celli]/V())*T)*diskDir_;
+        Usource[celli] += ((cellsV[celli]/V())*T)*diskDir;
     }
 
     if
@@ -128,7 +129,7 @@ void Foam::fv::actuationDiskSource::calcFroudeMethod
         writeCurrentTime(os);
 
         os  << Uref << tab << Cp << tab << Ct << tab << a << tab << T
-            << endl;
+            << tab << diskDir << endl;
     }
 }
 
@@ -221,12 +222,13 @@ void Foam::fv::actuationDiskSource::calcVariableScalingMethod
     const scalar CpStar = Cp*pow3(magUref/magUdisk);
 
     // Compute calibrated thrust/power (LSRMTK:Eq. 5)
-    const scalar T = 0.5*rhoRef*diskArea_*magSqr(Udisk & diskDir_)*CtStar;
-    const scalar P = 0.5*rhoRef*diskArea_*pow3(mag(Udisk & diskDir_))*CpStar;
+    const vector diskDir = this->diskDir();
+    const scalar T = 0.5*rhoRef*diskArea_*magSqr(Udisk & diskDir)*CtStar;
+    const scalar P = 0.5*rhoRef*diskArea_*pow3(mag(Udisk & diskDir))*CpStar;
 
     for (const label celli : cells_)
     {
-        Usource[celli] += (cellsV[celli]/totalV*T)*diskDir_;
+        Usource[celli] += (cellsV[celli]/totalV*T)*diskDir;
     }
 
     if
@@ -240,7 +242,7 @@ void Foam::fv::actuationDiskSource::calcVariableScalingMethod
 
         os  << Uref << tab << Cp << tab << Ct << tab
             << Udisk << tab << CpStar << tab << CtStar << tab << T << tab << P
-            << endl;
+            << tab << diskDir << endl;
     }
 }
 

src/fvOptions/sources/derived/radialActuationDiskSource/radialActuationDiskSourceTemplates.C

@@ -48,9 +48,8 @@ addRadialActuationDiskAxialInertialResistance
     scalarField Tr(cells.size());
 
     tensor E(Zero);
-    E.xx() = diskDir_.x();
-    E.yy() = diskDir_.y();
-    E.zz() = diskDir_.z();
+    const vector diskDir = this->diskDir();
+    E.diag(diskDir);
 
     const Field<vector> zoneCellCentres(mesh().cellCentres(), cells);
     const Field<scalar> zoneCellVolumes(mesh().cellVolumes(), cells);

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/Make/files

@@ -216,6 +216,7 @@ $(pointSources)/searchableSurfaceToPoint/searchableSurfaceToPoint.C
 $(pointSources)/sphereToPoint/sphereToPoint.C
 $(pointSources)/surfaceToPoint/surfaceToPoint.C
 $(pointSources)/zoneToPoint/zoneToPoint.C
+$(pointSources)/patchToPoint/patchToPoint.C
 
 faceZoneSources = topoSet/faceZoneSources
 $(faceZoneSources)/topoSetFaceZoneSource/topoSetFaceZoneSource.C