ENH: pointMesh: support for subsetMesh

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(),
+                pMesh.thisDb().time().timeName(),
-                baseMesh,
+                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);
+                const pointMesh& pMesh = pointMesh::New
-                PtrList<pointMesh> pMeshes(procMeshes.meshes().size());
+                (
-
+                    mesh,
-                forAll(pMeshes, proci)
+                    IOobject::READ_IF_PRESENT
-                {
+                );
-                    pMeshes.set
-                    (
-                        proci,
-                        new pointMesh(procMeshes.meshes()[proci])
-                    );
-                }
 
                 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/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/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/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
+        //- Construct from components
-        explicit pointPatch(const pointBoundaryMesh& bm)
+        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/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/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/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
+    forAll(pmUbf, patchi)
-    cmUbf.setSize(fvMesh_.boundary().size());
-
-    forAll(cmUbf, 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)
+            if (isA<fixedValuePointPatchField<Type>>(pfld))
-        {
+            {
-            Pout<< "Patch:" << fvMesh_.boundary()[patchi].patch().name()
+                cmUbf[polyPatchi] = cellMotionFvPatchField<Type>::typeName;
-                << " pointType:" << pmUbf.types()[patchi]
+            }
-                << " cellType:" << cmUbf[patchi] << endl;
+            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 /= magD;
-            {
+            d *= min(magD, mag(clipVelocity));
-                d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
-            }
-            else
-            {
-                d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
-            }
         }
     }
 
     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.writeEntryIfDifferent<vector>
-    os.writeEntry("wedgePlane", wedgePlane_);
+    (
-
+        "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.writeEntryIfDifferent<vector>
-    os.writeEntry("wedgePlane", wedgePlane_);
+    (
-
+        "projectDirection",
+        Zero,
+        projectDir_
+    );
+    os.writeEntryIfDifferent<label>("wedgePlane", -1, wedgePlane_);
     os.writeEntryIfDifferent<word>
     (
         "frozenPointsZone",

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

src/meshTools/topoSet/pointSources/patchToPoint/patchToPoint.C

@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "patchToPoint.H"
+#include "pointMesh.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(patchToPoint, 0);
+    addToRunTimeSelectionTable(topoSetSource, patchToPoint, word);
+    addToRunTimeSelectionTable(topoSetSource, patchToPoint, istream);
+    addToRunTimeSelectionTable(topoSetPointSource, patchToPoint, word);
+    addToRunTimeSelectionTable(topoSetPointSource, patchToPoint, istream);
+    addNamedToRunTimeSelectionTable
+    (
+        topoSetPointSource,
+        patchToPoint,
+        word,
+        patch
+    );
+    addNamedToRunTimeSelectionTable
+    (
+        topoSetPointSource,
+        patchToPoint,
+        istream,
+        patch
+    );
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::patchToPoint::usage_
+(
+    patchToPoint::typeName,
+    "\n    Usage: patchToPoint patch\n\n"
+    "    Select all points in the pointPatch."
+    " Note:accepts wildcards for patch.\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::patchToPoint::combine(topoSet& set, const bool add) const
+{
+    const pointMesh& pMesh = pointMesh::New(mesh_, IOobject::READ_IF_PRESENT);
+    const pointBoundaryMesh& pbm = pMesh.boundary();
+
+    labelList patchIDs
+    (
+        pbm.indices
+        (
+            selectedPatches_,
+            true        // useGroups
+        )
+    );
+
+
+    for (const label patchi : patchIDs)
+    {
+        const pointPatch& pp = pbm[patchi];
+
+        if (verbose_)
+        {
+            Info<< "    Found matching patch " << pp.name() << " with "
+                << returnReduce(pp.size(), sumOp<label>()) << " points" << endl;
+        }
+
+        for (const label pointi : pp.meshPoints())
+        {
+            addOrDelete(set, pointi, add);
+        }
+    }
+
+    if (patchIDs.empty())
+    {
+        WarningInFunction
+            << "Cannot find any patches matching "
+            << flatOutput(selectedPatches_) << nl
+            //<< "Valid names: " << flatOutput(mesh_.boundaryMesh().names())
+            << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::patchToPoint::patchToPoint
+(
+    const polyMesh& mesh,
+    const wordRe& patchName
+)
+:
+    topoSetPointSource(mesh),
+    selectedPatches_(one{}, patchName)
+{}
+
+
+Foam::patchToPoint::patchToPoint
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetPointSource(mesh),
+    selectedPatches_()
+{
+    // Look for 'patches' and 'patch', but accept 'name' as well
+    if (!dict.readIfPresent("patches", selectedPatches_))
+    {
+        selectedPatches_.resize(1);
+        selectedPatches_.front() =
+            dict.getCompat<wordRe>("patch", {{"name", 1806}});
+    }
+}
+
+
+Foam::patchToPoint::patchToPoint
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetPointSource(mesh),
+    selectedPatches_(one{}, wordRe(checkIs(is)))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::patchToPoint::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if (action == topoSetSource::ADD || action == topoSetSource::NEW)
+    {
+        if (verbose_)
+        {
+            Info<< "    Adding all points of patches: "
+                << flatOutput(selectedPatches_) << " ..." << endl;
+        }
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::SUBTRACT)
+    {
+        if (verbose_)
+        {
+            Info<< "    Removing all points of patches: "
+                << flatOutput(selectedPatches_) << " ..." << endl;
+        }
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/topoSet/pointSources/patchToPoint/patchToPoint.H

@@ -0,0 +1,174 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::patchToPoint
+
+Description
+    A \c topoSetPointSource to select points associated with given patch(es).
+
+    Operands:
+    \table
+      Operand   | Type     | Location
+      output    | pointSet | $FOAM_CASE/constant/polyMesh/sets/\<set\>
+    \endtable
+
+Usage
+    Minimal example by using \c system/topoSetDict.actions:
+    \verbatim
+    {
+        // Mandatory (inherited) entries
+        name        <name>;
+        type        pointSet;
+        action      <action>;
+
+        // Mandatory entries
+        source      patchToPoint;
+
+        // Conditional mandatory entries
+        // Select either of the below
+
+        // Option-1
+        patches
+        (
+            <patchName1>
+            <patchName2>
+            ...
+        );
+
+        // Option-2
+        patch    <patchName>;
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                         | Type | Req'd | Dflt
+      name       | Name of pointSet                    | word |  yes  | -
+      type       | Type name: pointSet                 | word |  yes  | -
+      action     | Action applied on points - see below | word |  yes  | -
+      source     | Source name: patchToPoint            | word |  yes  | -
+    \endtable
+
+    Options for the \c action entry:
+    \verbatim
+      new      | Create a new pointSet from selected points
+      add      | Add selected points into this pointSet
+      subtract | Remove selected points from this pointSet
+    \endverbatim
+
+    Options for the conditional mandatory entries:
+    \verbatim
+      Entry    | Description           | Type     | Req'd  | Dflt
+      patches  | Names of patches      | wordList | cond'l | -
+      patch    | Name of patch         | word     | cond'l | -
+    \endverbatim
+
+Note
+    The order of precedence among the conditional mandatory entries from the
+    highest to the lowest is \c patches, and \c patch.
+
+See also
+    - Foam::patchToFace
+
+SourceFiles
+    patchToPoint.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef patchToPoint_H
+#define patchToPoint_H
+
+#include "topoSetPointSource.H"
+#include "wordRes.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class patchToPoint Declaration
+\*---------------------------------------------------------------------------*/
+
+class patchToPoint
+:
+    public topoSetPointSource
+{
+    // Private Data
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Matcher for patches
+        wordRes selectedPatches_;
+
+
+    // Private Member Functions
+
+        void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("patchToPoint");
+
+
+    // Constructors
+
+        //- Construct from components
+        patchToPoint(const polyMesh& mesh, const wordRe& patchName);
+
+        //- Construct from dictionary
+        patchToPoint(const polyMesh& mesh, const dictionary& dict);
+
+        //- Construct from Istream
+        patchToPoint(const polyMesh& mesh, Istream& is);
+
+
+    //- Destructor
+    virtual ~patchToPoint() = default;
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet& set
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/parallel/reconstruct/reconstruct/processorMeshes.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.
@@ -31,6 +31,7 @@ License
 #include "IndirectList.H"
 #include "primitiveMesh.H"
 #include "OSspecific.H"
+#include "pointMesh.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -46,6 +47,9 @@ void Foam::processorMeshes::read()
 {
     // Make sure to clear (and hence unregister) any previously loaded meshes
     // and fields
+    pBoundaryProcAddressing_.free();
+    pMeshes_.free();
+
     boundaryProcAddressing_.free();
     cellProcAddressing_.free();
     faceProcAddressing_.free();
@@ -92,6 +96,39 @@ void Foam::processorMeshes::read()
         // boundaryProcAddressing (polyMesh)
         ioAddr.rename("boundaryProcAddressing");
         boundaryProcAddressing_.emplace_set(proci, ioAddr);
+
+
+        // pointMesh
+        // ~~~~~~~~~
+
+        pMeshes_.set
+        (
+            proci,
+            new pointMesh
+            (
+                meshes_[proci],
+                IOobject::READ_IF_PRESENT
+            )
+        );
+
+        pBoundaryProcAddressing_.set
+        (
+            proci,
+            autoPtr<labelIOList>::New
+            (
+                IOobject
+                (
+                    "boundaryProcAddressing",
+                    meshes_[proci].facesInstance(),
+                    polyMesh::meshSubDir/pointMesh::meshSubDir,
+                    pMeshes_[proci].thisDb(),
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                boundaryProcAddressing_[proci]
+            )
+        );
     }
 }
 
@@ -110,7 +147,9 @@ Foam::processorMeshes::processorMeshes
     pointProcAddressing_(databases.size()),
     faceProcAddressing_(databases.size()),
     cellProcAddressing_(databases.size()),
-    boundaryProcAddressing_(databases.size())
+    boundaryProcAddressing_(databases.size()),
+    pMeshes_(databases.size()),
+    pBoundaryProcAddressing_(databases.size())
 {
     read();
 }
@@ -250,6 +289,26 @@ void Foam::processorMeshes::removeFiles(const polyMesh& mesh)
     // boundaryProcAddressing
     io.rename("boundaryProcAddressing");
     fileHandler().rm(fileHandler().filePath(io.objectPath()));
+
+
+
+    // pointMesh
+    // ~~~~~~~~~
+
+    IOobject pointIO
+    (
+        "boundary",
+        mesh.facesInstance(),
+        polyMesh::meshSubDir/pointMesh::meshSubDir,
+        mesh.thisDb()
+    );
+
+    // pointMesh/boundary
+    fileHandler().rm(fileHandler().filePath(pointIO.objectPath()));
+
+    // boundaryProcAddressing
+    io.rename("boundaryProcAddressing");
+    fileHandler().rm(fileHandler().filePath(pointIO.objectPath()));
 }
 
 

src/parallel/reconstruct/reconstruct/processorMeshes.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2014 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2016,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -77,6 +77,15 @@ class processorMeshes
         PtrList<labelIOList> boundaryProcAddressing_;
 
 
+        // pointMesh
+
+            //- List of processor meshes
+            PtrList<pointMesh> pMeshes_;
+
+            //- List of processor boundary addressing lists for pointMeshes
+            PtrList<labelIOList> pBoundaryProcAddressing_;
+
+
     // Private Member Functions
 
         //- Read all meshes
@@ -140,6 +149,20 @@ public:
         }
 
 
+        // pointMesh
+
+            const PtrList<pointMesh>& pointMeshes() const noexcept
+            {
+                return pMeshes_;
+            }
+
+            const PtrList<labelIOList>& pointMeshBoundaryProcAddressing()
+            const noexcept
+            {
+                return pBoundaryProcAddressing_;
+            }
+
+
         //- Helper: remove all procAddressing files from mesh instance
         static void removeFiles(const polyMesh& mesh);
 };

tutorials/mesh/blockMesh/sphere7/0/pointDisplacement

@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2309                                  |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    arch        "LSB;label=32;scalar=64";
+    class       pointVectorField;
+    location    "0";
+    object      pointDisplacement;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 0 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    minX
+    {
+        type            surfaceDisplacement;
+        value           uniform (0 0 0);
+        velocity        (0.1 0.1 0.1);
+        geometry
+        {
+            box
+            {
+                type            triSurfaceMesh;
+                file            box222.obj;
+            }
+        }
+        projectMode     nearest;
+    }
+    maxX
+    {
+        $minX
+    }
+    minY
+    {
+        $minX
+    }
+    maxY
+    {
+        $minX
+    }
+    minZ
+    {
+        $minX
+    }
+    maxZ
+    {
+        $minX
+    }
+    boundaryEdges
+    {
+        type            edgeSlipDisplacement;
+        file            "box222.eMesh";
+        //file            "box222.extendedFeatureEdgeMesh";
+        //- Underrelax motion by limiting displacement velocity
+        velocity        (0.1 0.1 0.1);
+    }
+    boundaryPoints
+    {
+        type            pointAttraction;
+        file            "box222.eMesh";
+        //- Underrelax motion by limiting displacement velocity
+        velocity        (0.1 0.1 0.1);
+    }
+}
+
+
+// ************************************************************************* //

tutorials/mesh/blockMesh/sphere7/Allclean

@@ -3,8 +3,15 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
 #------------------------------------------------------------------------------
 
-cleanCase0
+cleanCase
 
-rm -rf constant
+rm -rf constant/extendedFeatureEdgeMesh
+
+#- From pointMesh generation
+rm -rf constant/pointMesh
+rm -f constant/triSurface/blockMesh.obj
+
+#- From surfaceFeatureExtract
+rm -f constant/triSurface/box222.eMesh
 
 #------------------------------------------------------------------------------

tutorials/mesh/blockMesh/sphere7/Allrun

@@ -3,6 +3,31 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
+#- Run blockMesh with projection
 runApplication blockMesh
 
+#- Generate pointMesh with additional feature patches
+#  (in constant/pointMesh/boundary). Extracted outside of mesh is not used.
+runApplication surfaceMeshExtract -featureAngle 45 \
+    constant/triSurface/blockMesh.obj
+
+
+#- TBD. move triSurfaces to resources
+#mkdir -p constant/triSurface
+#cp -f \
+#    "$FOAM_TUTORIALS"/resources/geometry/box222.obj \
+#    constant/triSurface/
+
+#- Extract features from surface. Writes .eMesh
+runApplication surfaceFeatureExtract
+
+##- Morph mesh to the surface (.obj file) and feature-edges (.eMesh file)
+#runApplication moveDynamicMesh
+
+runApplication decomposePar
+
+runParallel moveDynamicMesh
+
+runApplication reconstructPar
+
 #------------------------------------------------------------------------------

tutorials/mesh/blockMesh/sphere7/constant/dynamicMeshDict

@@ -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;
+    object      dynamicMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dynamicFvMesh   dynamicMotionSolverFvMesh;
+
+motionSolverLibs (fvMotionSolvers);
+
+motionSolver displacementSBRStress;
+diffusivity  uniform;
+
+// ************************************************************************* //

tutorials/mesh/blockMesh/sphere7/constant/triSurface/box222.obj

@@ -0,0 +1,79 @@
+# Generated by Visualization Toolkit
+v -1 -1 -1
+v -1 -1 1
+v -1 1 -1
+v -1 1 1
+v 1 -1 -1
+v 1 -1 1
+v 1 1 -1
+v 1 1 1
+v -1 -1 -1
+v -1 -1 1
+v 1 -1 -1
+v 1 -1 1
+v -1 1 -1
+v -1 1 1
+v 1 1 -1
+v 1 1 1
+v -1 -1 -1
+v 1 -1 -1
+v -1 1 -1
+v 1 1 -1
+v -1 -1 1
+v 1 -1 1
+v -1 1 1
+v 1 1 1
+vn -1 0 0
+vn -1 0 0
+vn -1 0 0
+vn -1 0 0
+vn 1 0 0
+vn 1 0 0
+vn 1 0 0
+vn 1 0 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 1 0
+vn 0 1 0
+vn 0 1 0
+vn 0 1 0
+vn 0 0 -1
+vn 0 0 -1
+vn 0 0 -1
+vn 0 0 -1
+vn 0 0 1
+vn 0 0 1
+vn 0 0 1
+vn 0 0 1
+vt -0.5 -0.5
+vt 1.5 -0.5
+vt -0.5 1.5
+vt 1.5 1.5
+vt 0.5 -0.5
+vt -1.5 -0.5
+vt 0.5 1.5
+vt -1.5 1.5
+vt 0.5 0.5
+vt 0.5 -1.5
+vt -1.5 0.5
+vt -1.5 -1.5
+vt -0.5 0.5
+vt -0.5 -1.5
+vt 1.5 0.5
+vt 1.5 -1.5
+vt 0.5 -0.5
+vt -1.5 -0.5
+vt 0.5 1.5
+vt -1.5 1.5
+vt -0.5 -0.5
+vt 1.5 -0.5
+vt -0.5 1.5
+vt 1.5 1.5
+f 1/1/1 2/2/2 4/4/4 3/3/3
+f 5/5/5 7/7/7 8/8/8 6/6/6
+f 9/9/9 11/11/11 12/12/12 10/10/10
+f 13/13/13 14/14/14 16/16/16 15/15/15
+f 17/17/17 19/19/19 20/20/20 18/18/18
+f 21/21/21 22/22/22 24/24/24 23/23/23

tutorials/mesh/blockMesh/sphere7/constant/triSurface/box422.obj

@@ -0,0 +1,79 @@
+# Generated by Visualization Toolkit
+v -2 -1 -1
+v -2 -1 1
+v -2 1 -1
+v -2 1 1
+v 2 -1 -1
+v 2 -1 1
+v 2 1 -1
+v 2 1 1
+v -2 -1 -1
+v -2 -1 1
+v 2 -1 -1
+v 2 -1 1
+v -2 1 -1
+v -2 1 1
+v 2 1 -1
+v 2 1 1
+v -2 -1 -1
+v 2 -1 -1
+v -2 1 -1
+v 2 1 -1
+v -2 -1 1
+v 2 -1 1
+v -2 1 1
+v 2 1 1
+vn -1 0 0
+vn -1 0 0
+vn -1 0 0
+vn -1 0 0
+vn 1 0 0
+vn 1 0 0
+vn 1 0 0
+vn 1 0 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 -1 0
+vn 0 1 0
+vn 0 1 0
+vn 0 1 0
+vn 0 1 0
+vn 0 0 -1
+vn 0 0 -1
+vn 0 0 -1
+vn 0 0 -1
+vn 0 0 1
+vn 0 0 1
+vn 0 0 1
+vn 0 0 1
+vt -0.5 -0.5
+vt 1.5 -0.5
+vt -0.5 1.5
+vt 1.5 1.5
+vt 0.5 -0.5
+vt -1.5 -0.5
+vt 0.5 1.5
+vt -1.5 1.5
+vt 1.5 0.5
+vt 1.5 -1.5
+vt -2.5 0.5
+vt -2.5 -1.5
+vt -1.5 0.5
+vt -1.5 -1.5
+vt 2.5 0.5
+vt 2.5 -1.5
+vt 1.5 -0.5
+vt -2.5 -0.5
+vt 1.5 1.5
+vt -2.5 1.5
+vt -1.5 -0.5
+vt 2.5 -0.5
+vt -1.5 1.5
+vt 2.5 1.5
+f 1/1/1 2/2/2 4/4/4 3/3/3
+f 5/5/5 7/7/7 8/8/8 6/6/6
+f 9/9/9 11/11/11 12/12/12 10/10/10
+f 13/13/13 14/14/14 16/16/16 15/15/15
+f 17/17/17 19/19/19 20/20/20 18/18/18
+f 21/21/21 22/22/22 24/24/24 23/23/23

tutorials/mesh/blockMesh/sphere7/system/blockMeshDict

@@ -127,16 +127,56 @@ faces
 
 boundary
 (
-    walls
+    minX
     {
-        type wall;
+        type patch;
         faces
         (
             (1 0)  // block 1: x-min
+        );
+    }
+
+    maxX
+    {
+        type patch;
+        faces
+        (
             (2 1)  // block 2: x-max
+        );
+    }
+
+    minY
+    {
+        type patch;
+        faces
+        (
             (3 2)  // block 3: y-min
+        );
+    }
+
+    maxY
+    {
+        type patch;
+        faces
+        (
             (4 3)  // block 4: y-max
+        );
+    }
+
+    minZ
+    {
+        type patch;
+        faces
+        (
             (5 4)  // block 5: z-min
+        );
+    }
+
+    maxZ
+    {
+        type patch;
+        faces
+        (
             (6 5)  // block 6: z-max
         );
     }

tutorials/mesh/blockMesh/sphere7/system/controlDict

@@ -14,6 +14,12 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+DebugSwitches
+{
+    pointBoundaryMesh   1;
+    fvMeshSubset        1;
+}
+
 application     blockMesh;
 
 startFrom       startTime;
@@ -22,9 +28,9 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0;
+endTime         10;
 
-deltaT          0;
+deltaT          1;
 
 writeControl    timeStep;
 

tutorials/mesh/blockMesh/sphere7/system/fvSchemes

@@ -18,19 +18,27 @@ ddtSchemes
 {}
 
 gradSchemes
-{}
+{
+    default         Gauss linear;
+}
 
 divSchemes
 {}
 
 laplacianSchemes
-{}
+{
+    default         Gauss linear corrected;
+}
 
 interpolationSchemes
-{}
+{
+    default         linear;
+}
 
 snGradSchemes
-{}
+{
+    default         corrected;
+}
 
 
 // ************************************************************************* //

tutorials/mesh/blockMesh/sphere7/system/fvSolution

@@ -14,5 +14,20 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+solvers
+{
+    "cellDisplacement.*"
+    {
+        solver          GAMG;
+        tolerance       1e-08;
+        relTol          0;
+        smoother        GaussSeidel;
+    }
+}
+
+
+PIMPLE
+{}
+
 
 // ************************************************************************* //

tutorials/mesh/blockMesh/sphere7/system/surfaceFeatureExtractDict

@@ -0,0 +1,33 @@
+/*--------------------------------*- 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      surfaceFeatureExtractDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+box222.obj
+{
+    // How to obtain raw features (extractFromFile || extractFromSurface)
+    extractionMethod    extractFromSurface;
+
+    // Mark edges whose adjacent surface normals are at an angle less
+    // than includedAngle as features
+    // - 0  : selects no edges
+    // - 180: selects all edges
+    includedAngle       135;
+
+    // Write features to obj format for postprocessing
+    writeObj            yes;
+}
+
+
+// ************************************************************************* //