Revert "Merge branch 'primitiveMeshOptimization' into 'develop'"

This reverts merge request !596

applications/test/boundBox2/Test-boundBox2.C

@@ -143,7 +143,7 @@ int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
 
-    treeBoundBox bb(cube(0, 1));
+    treeBoundBox bb(zero_one{});
     treeBoundBox sub(cube(0.1, 0.8));
 
     Info<< nl

applications/test/polyMeshGeom-speed1/Make/files

@@ -0,0 +1,3 @@
+Test-polyMeshGeom-speed1.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-polyMeshGeom-speed1

applications/test/polyMeshGeom-speed1/Make/options

@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/mesh/blockMesh/lnInclude
+
+EXE_LIBS = \
+    -lblockMesh

applications/test/polyMeshGeom-speed1/Test-polyMeshGeom-speed1.C

@@ -0,0 +1,804 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+Description
+    Simple timing tests for some polyMesh primitives
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "clockTime.H"
+#include "Time.H"
+#include "PDRblock.H"
+#include "polyMesh.H"
+#include "ListOps.H"
+
+using namespace Foam;
+
+void printAlloc(const polyMesh& mesh)
+{
+    Info<< "memory"
+        << " hasCellPoints:" << mesh.hasCellPoints()
+        << " hasPointCells:" << mesh.hasPointCells() << endl;
+}
+
+
+void printInfo(const polyMesh& mesh)
+{
+    Info<< "polyMesh"
+        << " nPoints:" << mesh.nPoints()
+        << " nInternalFaces:" << mesh.nInternalFaces()
+        << " nFaces:" << mesh.nFaces()
+        << " nCells:" << mesh.nCells() << endl;
+}
+
+
+// How point cells are calculated in OpenFOAM-v2212 and earlier
+autoPtr<labelListList> pointCells_2212(const polyMesh& mesh)
+{
+    const cellList& cf = mesh.cells();
+
+    // Count number of cells per point
+
+    labelList npc(mesh.nPoints(), Zero);
+
+    forAll(cf, celli)
+    {
+        const labelList curPoints = cf[celli].labels(mesh.faces());
+
+        for (const label pointi : curPoints)
+        {
+            ++npc[pointi];
+        }
+    }
+
+
+    // Size and fill cells per point
+
+    auto pcPtr_ = autoPtr<labelListList>::New(npc.size());
+    labelListList& pointCellAddr = *pcPtr_;
+
+    forAll(pointCellAddr, pointi)
+    {
+        pointCellAddr[pointi].setSize(npc[pointi]);
+        npc[pointi] = 0;
+    }
+
+
+    forAll(cf, celli)
+    {
+        const labelList curPoints = cf[celli].labels(mesh.faces());
+
+        for (const label pointi : curPoints)
+        {
+            pointCellAddr[pointi][npc[pointi]++] = celli;
+        }
+    }
+
+    return pcPtr_;
+}
+
+
+// Line cell::labels but with persistent storage
+void cell_labels
+(
+    const cell& cFaces,
+    const faceUList& meshFaces,
+    DynamicList<label>& pointLabels
+)
+{
+    // const labelList& cFaces = *this;
+
+    label nVerts = 0;
+    for (const label facei : cFaces)
+    {
+        nVerts += meshFaces[facei].size();
+    }
+
+    // pointLabels.clear();
+    pointLabels.expandStorage();
+
+    // The first face has no duplicates, can copy in values
+    const labelList& firstFace = meshFaces[cFaces[0]];
+
+    std::copy(firstFace.cbegin(), firstFace.cend(), pointLabels.begin());
+
+    // Now already contains some vertices
+    nVerts = firstFace.size();
+
+    // For the rest of the faces. For each vertex, check if the point is
+    // already inserted (up to nVerts, which now carries the number of real
+    // points. If not, add it at the end of the list.
+
+    for (label facei = 1; facei < cFaces.size(); ++facei)
+    {
+        for (const label curPoint : meshFaces[cFaces[facei]])
+        {
+            bool pointFound = false;
+
+            for (label checki = 0; checki < nVerts; ++checki)
+            {
+                if (curPoint == pointLabels[checki])
+                {
+                    pointFound = true;
+                    break;
+                }
+            }
+
+            if (!pointFound)
+            {
+                pointLabels[nVerts] = curPoint;
+                ++nVerts;
+            }
+        }
+    }
+
+    pointLabels.resize(nVerts);
+}
+
+
+
+// Like OpenFOAM-v2212, but with cell::labels unrolled to avoid allocations
+autoPtr<labelListList> pointCells_2212mod(const polyMesh& mesh)
+{
+    const cellList& cf = mesh.cells();
+
+    // Vertex labels for the current cell
+    DynamicList<label> vertices(256);
+
+    // Count number of cells per point
+
+    labelList npc(mesh.nPoints(), Zero);
+
+    for (const cell& c : cf)
+    {
+        cell_labels(c, mesh.faces(), vertices);
+
+        for (const label pointi : vertices)
+        {
+            ++npc[pointi];
+        }
+    }
+
+
+    // Size and fill cells per point
+
+    auto pcPtr_ = autoPtr<labelListList>::New(npc.size());
+    labelListList& pointCellAddr = *pcPtr_;
+
+    forAll(pointCellAddr, pointi)
+    {
+        pointCellAddr[pointi].resize(npc[pointi]);
+        npc[pointi] = 0;
+    }
+
+
+    forAll(cf, celli)
+    {
+        cell_labels(cf[celli], mesh.faces(), vertices);
+
+        for (const label pointi : vertices)
+        {
+            pointCellAddr[pointi][npc[pointi]++] = celli;
+        }
+    }
+
+    return pcPtr_;
+}
+
+
+// How cells points are calculated in OpenFOAM-v2212 and earlier
+autoPtr<labelListList> cellPoints_2212(const polyMesh& mesh)
+{
+    autoPtr<labelListList> pointCells = pointCells_2212(mesh);
+
+    auto cpPtr_ = autoPtr<labelListList>::New(mesh.nCells());
+
+    invertManyToMany(mesh.nCells(), pointCells(), *cpPtr_);
+
+    return cpPtr_;
+}
+
+
+// Calculate with bitSet tracking and avoid cells::labels
+autoPtr<labelListList> pointCells_bitSet(const polyMesh& mesh)
+{
+    // Calculate point-cell topology
+
+    const cellList& cellLst = mesh.cells();
+    const faceList& faceLst = mesh.faces();
+
+    // For tracking (only use each point id once)
+    bitSet usedPoints(mesh.nPoints());
+
+    // Vertex labels for the current cell
+    DynamicList<label> vertices(256);
+
+    const label loopLen = mesh.nCells();
+
+    // Step 1: count number of cells per point
+
+    labelList pointCount(mesh.nPoints(), Zero);
+
+    for (label celli = 0; celli < loopLen; ++celli)
+    {
+        // Clear any previous contents
+        usedPoints.unset(vertices);
+        vertices.clear();
+
+        for (const label facei : cellLst[celli])
+        {
+            for (const label pointi : faceLst[facei])
+            {
+                // Only once for each point id
+                if (usedPoints.set(pointi))
+                {
+                    vertices.push_back(pointi);
+                    ++pointCount[pointi];
+                }
+            }
+        }
+    }
+
+
+    // Step 2: set sizing, reset counters
+
+    auto pcPtr_ = autoPtr<labelListList>::New(mesh.nPoints());
+    auto& pointCellAddr = *pcPtr_;
+
+    forAll(pointCellAddr, pointi)
+    {
+        pointCellAddr[pointi].resize_nocopy(pointCount[pointi]);
+        pointCount[pointi] = 0;
+    }
+
+
+    // Step 3: fill in values. Logic as per step 1
+    for (label celli = 0; celli < loopLen; ++celli)
+    {
+        // Clear any previous contents
+        usedPoints.unset(vertices);
+        vertices.clear();
+
+        for (const label facei : cellLst[celli])
+        {
+            for (const label pointi : faceLst[facei])
+            {
+                // Only once for each point id
+                if (usedPoints.set(pointi))
+                {
+                    vertices.push_back(pointi);
+                    pointCellAddr[pointi][pointCount[pointi]++] = celli;
+                }
+            }
+        }
+    }
+
+    return pcPtr_;
+}
+
+
+// Calculate with bitSet tracking and avoid cells::labels
+autoPtr<labelListList> cellPoints_bitSet(const polyMesh& mesh)
+{
+    // Calculate cell-point topology
+
+    auto cpPtr_ = autoPtr<labelListList>::New(mesh.nCells());
+    auto& cellPointAddr = *cpPtr_;
+
+    const cellList& cellLst = mesh.cells();
+    const faceList& faceLst = mesh.faces();
+
+    // For tracking (only use each point id once)
+    bitSet usedPoints(mesh.nPoints());
+
+    // Vertex labels for the current cell
+    DynamicList<label> vertices(256);
+
+    const label loopLen = mesh.nCells();
+
+    for (label celli = 0; celli < loopLen; ++celli)
+    {
+        // Clear any previous contents
+        usedPoints.unset(vertices);
+        vertices.clear();
+
+        for (const label facei : cellLst[celli])
+        {
+            for (const label pointi : faceLst[facei])
+            {
+                // Only once for each point id
+                if (usedPoints.set(pointi))
+                {
+                    vertices.push_back(pointi);
+                }
+            }
+        }
+
+        cellPointAddr[celli] = vertices;  // unsorted
+    }
+
+    return cpPtr_;
+}
+
+
+// Calculate with linear lookup and avoid cells::labels
+autoPtr<labelListList> pointCells_linear(const polyMesh& mesh)
+{
+    // Calculate point-cell topology
+
+    const cellList& cellLst = mesh.cells();
+    const faceList& faceLst = mesh.faces();
+
+    // Vertex labels for the current cell
+    DynamicList<label> vertices(256);
+
+    const label loopLen = mesh.nCells();
+
+    // Step 1: count number of cells per point
+
+    labelList pointCount(mesh.nPoints(), Zero);
+
+    for (label celli = 0; celli < loopLen; ++celli)
+    {
+        // Clear any previous contents
+        vertices.clear();
+
+        for (const label facei : cellLst[celli])
+        {
+            for (const label pointi : faceLst[facei])
+            {
+                // Only once for each point id
+                if (!vertices.contains(pointi))
+                {
+                    vertices.push_back(pointi);
+                    ++pointCount[pointi];
+                }
+            }
+        }
+    }
+
+
+    // Step 2: set sizing, reset counters
+
+    auto pcPtr_ = autoPtr<labelListList>::New(mesh.nPoints());
+    auto& pointCellAddr = *pcPtr_;
+
+    forAll(pointCellAddr, pointi)
+    {
+        pointCellAddr[pointi].resize_nocopy(pointCount[pointi]);
+        pointCount[pointi] = 0;
+    }
+
+
+    // Step 3: fill in values. Logic as per step 1
+    for (label celli = 0; celli < loopLen; ++celli)
+    {
+        // Clear any previous contents
+        vertices.clear();
+
+        for (const label facei : cellLst[celli])
+        {
+            for (const label pointi : faceLst[facei])
+            {
+                // Only once for each point id
+                if (!vertices.contains(pointi))
+                {
+                    vertices.push_back(pointi);
+                    pointCellAddr[pointi][pointCount[pointi]++] = celli;
+                }
+            }
+        }
+    }
+
+    return pcPtr_;
+}
+
+
+// Calculate with linear lookup and avoid cells::labels
+autoPtr<labelListList> cellPoints_linear(const polyMesh& mesh)
+{
+    // Calculate cell-point topology
+
+    auto cpPtr_ = autoPtr<labelListList>::New(mesh.nCells());
+    auto& cellPointAddr = *cpPtr_;
+
+    const cellList& cellLst = mesh.cells();
+    const faceList& faceLst = mesh.faces();
+
+    // Vertex labels for the current cell
+    DynamicList<label> vertices(256);
+
+    const label loopLen = mesh.nCells();
+
+    for (label celli = 0; celli < loopLen; ++celli)
+    {
+        // Clear any previous contents
+        vertices.clear();
+
+        for (const label facei : cellLst[celli])
+        {
+            for (const label pointi : faceLst[facei])
+            {
+                // Only once for each point id
+                if (!vertices.contains(pointi))
+                {
+                    vertices.push_back(pointi);
+                }
+            }
+        }
+
+        cellPointAddr[celli] = vertices;  // unsorted
+    }
+
+    return cpPtr_;
+}
+
+
+// Calculate point-cell from point-face information
+autoPtr<labelListList> pointCells_faces(const polyMesh& mesh)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+    const labelListList& pFaces = mesh.pointFaces();
+
+    const label loopLen = mesh.nPoints();
+
+    auto pcPtr_ = autoPtr<labelListList>::New(mesh.nPoints());
+    auto& pointCellAddr = *pcPtr_;
+
+    DynamicList<label> storage(256);
+
+    for (label pointi = 0; pointi < loopLen; ++pointi)
+    {
+        // Clear any previous contents
+        storage.clear();
+
+        for (const label facei : pFaces[pointi])
+        {
+            // Owner cell
+            storage.push_back(own[facei]);
+
+            // Neighbour cell
+            if (facei < mesh.nInternalFaces())
+            {
+                storage.push_back(nei[facei]);
+            }
+        }
+
+        // Sort + unique to eliminate duplicates
+        std::sort(storage.begin(), storage.end());
+        auto last = std::unique(storage.begin(), storage.end());
+        storage.resize(label(last - storage.begin()));
+
+        pointCellAddr[pointi] = storage;
+    }
+
+    return pcPtr_;
+}
+
+// Calculate point-cell from point-face information
+autoPtr<labelListList> pointCells_bitSet_faces(const polyMesh& mesh)
+{
+    const labelList& own = mesh.faceOwner();
+    const labelList& nei = mesh.faceNeighbour();
+    const labelListList& pFaces = mesh.pointFaces();
+
+    const label loopLen = mesh.nPoints();
+
+    auto pcPtr_ = autoPtr<labelListList>::New(mesh.nPoints());
+    auto& pointCellAddr = *pcPtr_;
+
+    // For tracking (only use each cell id once)
+    bitSet usedCells(mesh.nCells());
+
+    DynamicList<label> storage(256);
+
+    for (label pointi = 0; pointi < loopLen; ++pointi)
+    {
+        // Clear any previous contents
+        usedCells.unset(storage);
+        storage.clear();
+
+        for (const label facei : pFaces[pointi])
+        {
+            // Owner cell - only once
+            if (usedCells.set(own[facei]))
+            {
+                storage.push_back(own[facei]);
+            }
+
+            // Neighbour cell
+            if (facei < mesh.nInternalFaces() && usedCells.set(nei[facei]))
+            {
+                storage.push_back(nei[facei]);
+            }
+        }
+
+        pointCellAddr[pointi] = storage;
+    }
+
+    return pcPtr_;
+}
+
+
+// Calculate point-cell from cell-point information
+autoPtr<labelListList> pointCells_bitSet_alon(const polyMesh& mesh)
+{
+    autoPtr<labelListList> cellPoints = cellPoints_bitSet(mesh);
+
+    auto pcPtr_ = autoPtr<labelListList>::New(mesh.nPoints());
+
+    invertManyToMany(mesh.nPoints(), cellPoints(), *pcPtr_);
+
+    return pcPtr_;
+}
+
+
+// Eliminate duplicates with sort+unique
+autoPtr<labelListList> cellPoints_sorted(const polyMesh& mesh)
+{
+    // Calculate cell-point topology
+
+    auto cpPtr_ = autoPtr<labelListList>::New(mesh.nCells());
+    auto& cellPointAddr = *cpPtr_;
+
+    const cellList& cellLst = mesh.cells();
+    const faceList& faceLst = mesh.faces();
+
+    // Vertex labels for the current cell
+    DynamicList<label> vertices(256);
+
+    const label loopLen = mesh.nCells();
+
+    for (label celli = 0; celli < loopLen; ++celli)
+    {
+        // Clear any previous contents
+        vertices.clear();
+
+        for (const label facei : cellLst[celli])
+        {
+            for (const label pointi : faceLst[facei])
+            {
+                vertices.push_back(pointi);
+            }
+        }
+
+        // Sort + unique to eliminate duplicates
+        std::sort(vertices.begin(), vertices.end());
+        auto last = std::unique(vertices.begin(), vertices.end());
+        vertices.resize(label(last - vertices.begin()));
+
+        cellPointAddr[celli] = vertices;
+    }
+
+    return cpPtr_;
+}
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    argList::noBanner();
+    argList::noParallel();
+    argList::noFunctionObjects();
+    argList::addOption("nCells", "number", "The number of cells");
+
+    #include "setRootCase.H"
+
+    const scalar cellCount(args.getOrDefault<scalar>("nCells", 1000));
+
+    const label nDivs(::round(::cbrt(cellCount)));
+
+    PDRblock blkMesh(boundBox(zero_one{}), labelVector::uniform(nDivs));
+
+    autoPtr<Time> dummyTimePtr(Time::New());
+
+    Info<< "Requested " << cellCount
+        << " cells, blockMesh with " << blkMesh.nCells() << " cells" << nl;
+
+    autoPtr<polyMesh> meshPtr = blkMesh.innerMesh
+    (
+        IOobject
+        (
+            "Testing",
+            dummyTimePtr->system(),
+            *dummyTimePtr,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            IOobject::NO_REGISTER
+        )
+    );
+
+    auto& mesh = meshPtr();
+
+    printInfo(mesh);
+    printAlloc(mesh);
+
+    clockTime timing;
+
+    // pointCells
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) mesh.pointCells();
+        Info<< "pointCells (builtin): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    // cellPoints
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) mesh.cellPoints();
+        Info<< "cellPoints (builtin): " << timing.elapsedTime() << " s" << nl;
+    }
+    Info<< nl;
+
+
+    // pointCells
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) pointCells_2212(mesh);
+        Info<< "pointCells (2212): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) pointCells_2212mod(mesh);
+        Info<< "pointCells (2212mod): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) pointCells_bitSet(mesh);
+        Info<< "pointCells (bitSet): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) pointCells_linear(mesh);
+        Info<< "pointCells (linear): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) pointCells_faces(mesh);
+        Info<< "pointCells (faces): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) pointCells_bitSet_faces(mesh);
+        Info<< "pointCells (bitSet faces): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) pointCells_bitSet_alon(mesh);
+        Info<< "pointCells (bitSet alon): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    // cellPoints
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) cellPoints_2212(mesh);
+        Info<< "cellPoints (2212): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) cellPoints_bitSet(mesh);
+        Info<< "cellPoints (bitSet): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) cellPoints_linear(mesh);
+        Info<< "cellPoints (linear): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) cellPoints_sorted(mesh);
+        Info<< "cellPoints (sorted): " << timing.elapsedTime() << " s" << nl;
+    }
+
+
+    // With precalculated values
+    {
+        mesh.clearOut();
+        const auto& cp = mesh.cellPoints();
+        timing.resetTime();
+
+        auto pcPtr_ = autoPtr<labelListList>::New(mesh.nPoints());
+
+        invertManyToMany(mesh.nPoints(), cp, *pcPtr_);
+
+        Info<< "pointCells (from cached cellPoints): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    // With precalculated values
+    {
+        mesh.clearOut();
+        (void)mesh.pointFaces();
+        timing.resetTime();
+
+        (void) pointCells_bitSet_faces(mesh);
+
+        Info<< "pointCells (bitSet from cached pointFaces): " << timing.elapsedTime() << " s" << nl;
+    }
+
+
+    // With precalculated values
+    {
+        mesh.clearOut();
+        const auto& pc = mesh.pointCells();
+        timing.resetTime();
+
+        auto cpPtr_ = autoPtr<labelListList>::New(mesh.nCells());
+
+        invertManyToMany(mesh.nCells(), pc, *cpPtr_);
+
+        Info<< "cellPoints (from cached pointCells): " << timing.elapsedTime() << " s" << nl;
+    }
+
+
+    // Re-measure timings
+    Info<< nl;
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) mesh.pointCells();
+        Info<< "pointCells (builtin): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    {
+        mesh.clearOut();
+        timing.resetTime();
+        (void) mesh.cellPoints();
+        Info<< "cellPoints (builtin): " << timing.elapsedTime() << " s" << nl;
+    }
+
+    Info<< "\nEnd\n" << nl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/triangleIntersection/Test-triangleIntersection.C

@@ -88,8 +88,7 @@ int main(int argc, char *argv[])
     // Info<<"tree-bb faces: " << treeBoundBox::faces << nl
     //     <<"tree-bb edges: " << treeBoundBox::edges << endl;
 
-    treeBoundBox bb;
-    bb = cube(0, 1);
+    treeBoundBox bb(zero_one{});
 
     triPoints tri;
     tri.a() = point(-0.1, 0.5, 0.5);

applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -76,9 +76,9 @@ Description
     writer.beginCellData(4);
     writer.writeProcIDs();
     {
-        // Use primitive patch order
         Field<scalar> fld
         (
+            // Use primitive patch order
             faMeshTools::flattenEdgeField(aMesh.magLe(), true)
         );
         writer.write("magLe", fld);
@@ -93,11 +93,17 @@ Description
 }
 
 {
+    const Field<vector> edgeCentres
+    (
+        // Use primitive patch order
+        faMeshTools::flattenEdgeField(aMesh.edgeCentres(), true)
+    );
+
     // finiteArea - edgeCentres
     // (no other convenient way to display vectors on the edges)
     vtk::lineWriter writer
     (
-        aMesh.edgeCentres(),
+        edgeCentres,
         edgeList::null(),
         // vtk::formatType::INLINE_ASCII,
         fileName
@@ -109,19 +115,20 @@ Description
     writer.writeGeometry();
 
     // PointData
-    writer.beginPointData(4);
+    writer.beginPointData(3);
+    writer.writeProcIDs();  // Unfortunately cannot threshold on points
     {
-        // Use primitive patch order
         Field<vector> fld
         (
+            // Use primitive patch order
             faMeshTools::flattenEdgeField(aMesh.Le(), true)
         );
         writer.write("Le", fld);
     }
     {
-        // Use primitive patch order
         Field<vector> fld
         (
+            // Use primitive patch order
             faMeshTools::flattenEdgeField(aMesh.edgeAreaNormals(), true)
         );
         writer.write("normal", fld);

bin/mpirunDebug

@@ -7,7 +7,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2015 OpenFOAM Foundation
-#     Copyright (C) 2017-2021 OpenCFD Ltd.
+#     Copyright (C) 2017-2023 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -45,6 +45,7 @@ options:
   -local        Same as -spawn=1
   -remote       Same as -spawn=2
   -clean        Remove old processor*.{log,sh} files, mpirun.schema etc
+  -decompose-dict=<file>   Specific decomposeParDict name
   -help         Print the usage
 
 Invoke mpirun with separate per-processor log files or running in
@@ -111,7 +112,7 @@ Linux)
 esac
 
 unset appName appArgs nProcs
-unset method spawn optClean
+unset method spawn optClean optValue
 optConfirm=true
 
 decompDict="system/decomposeParDict"
@@ -126,7 +127,7 @@ do
     then
         knownOption=true  # Assume success
         case "$1" in
-        '') ;;  # ignore junk
+        ('') ;;  # Ignore junk
 
         -clean)     optClean=true ;;
         -yes)       unset optConfirm ;;
@@ -173,14 +174,25 @@ do
             shift
             ;;
 
+        (-decompose-dict=*)
+            optValue="${1#*=}"
+            case "$optValue" in
+            ('' | none | false) ;;  ## Ignore
+            (*)
+                decompDict="$optValue"
+                appArgs="${appArgs}${appArgs:+ }-decomposeParDict '$decompDict'"
+                ;;
+            esac
+            ;;
+
         -decomposeParDict)
             # Grab values and add to args immediately
             decompDict="$2"
-            appArgs="${appArgs}${appArgs:+ }$1 '$2'"
             shift
+            appArgs="${appArgs}${appArgs:+ }-decomposeParDict '$decompDict'"
             ;;
 
-        *)
+        (*)
             knownOption=false  # Fallthrough to regular processing
             ;;
         esac
@@ -192,23 +204,24 @@ do
         fi
     fi
 
+    # Processing application arguments
     case "$1" in
-    -help* | --help*) usage ;;
-    '') ;;  # ignore junk
+    (-help* | --help*) usage ;;
+    ('') ;;  ## Ignore junk
 
-    -np)
+    (-np)
         nProcs="$2"
         shift
         ;;
 
-    -decomposeParDict)
+    (-decomposeParDict)
         # Grab values and add to args immediately
         decompDict="$2"
-        appArgs="${appArgs}${appArgs:+ }$1 '$2'"
+        appArgs="${appArgs}${appArgs:+ }-decomposeParDict '$decompDict'"
         shift
         ;;
 
-    *)
+    (*)
         if [ -z "$appName" ]
         then
             appName="$1"

bin/tools/RunFunctions

@@ -6,7 +6,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2015-2022 OpenCFD Ltd.
+#     Copyright (C) 2015-2023 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -180,6 +180,16 @@ getNumberOfProcessors()
 {
     local dict="${1:-system/decomposeParDict}"
 
+    case "$dict" in
+    (system/*)  # Already qualified
+        ;;
+    (*)
+        # If it does not exist, assume it refers to location in system/
+        [ -f "$dict" ] || dict="system/$dict"
+        ;;
+    esac
+
+
     # Re-use positional parameters for automatic whitespace elimination
     set -- $(foamDictionary -entry numberOfSubdomains -value "$dict" 2>/dev/null)
 
@@ -223,7 +233,7 @@ getApplication()
 #
 runApplication()
 {
-    local appName appRun logFile logMode
+    local appName appRun optValue logFile logMode
 
     # Any additional parsed arguments (eg, decomposeParDict)
     local appArgs
@@ -232,25 +242,39 @@ runApplication()
     while [ "$#" -gt 0 ] && [ -z "$appRun" ]
     do
         case "$1" in
-            -a | -append)
-                logMode=append
-                ;;
-            -o | -overwrite)
-                logMode=overwrite
-                ;;
-            -s | -suffix)
-                logFile=".$2"
-                shift
-                ;;
-            -decomposeParDict)
-                appArgs="$appArgs $1 $2"
-                shift
-                ;;
-            '')
-                ;;
-            *)
-                appRun="$1"
-                ;;
+        ('') ;;  # Ignore junk
+
+        (-a | -append)
+            logMode=append
+            ;;
+        (-o | -overwrite)
+            logMode=overwrite
+            ;;
+        (-s | -suffix)
+            logFile=".$2"
+            shift
+            ;;
+
+        (-decompose-dict=*)
+            optValue="${1#*=}"
+            case "$optValue" in
+            ('' | none | false) ;;  ## Ignore
+            (*) appArgs="$appArgs -decomposeParDict $optValue" ;;
+            esac
+            ;;
+
+        (-decomposeParDict)
+            optValue="$2"
+            shift
+            case "$optValue" in
+            ('' | none | false) ;; ## Ignore
+            (*) appArgs="$appArgs -decomposeParDict $optValue" ;;
+            esac
+            ;;
+
+        (*)
+            appRun="$1"
+            ;;
         esac
         shift
     done
@@ -280,7 +304,7 @@ runApplication()
 #
 runParallel()
 {
-    local appName appRun logFile logMode nProcs
+    local appName appRun optValue logFile logMode nProcs
 
     # Any additional parsed arguments (eg, decomposeParDict)
     local appArgs="-parallel"
@@ -295,30 +319,46 @@ runParallel()
     while [ "$#" -gt 0 ] && [ -z "$appRun" ]
     do
         case "$1" in
-            -a | -append)
-                logMode=append
+        ('') ;;  # Ignore junk
+
+        (-a | -append)
+            logMode=append
+            ;;
+        (-o | -overwrite)
+            logMode=overwrite
+            ;;
+        (-s | -suffix)
+            logFile=".$2"
+            shift
+            ;;
+        (-n | -np)
+            nProcs="$2"
+            shift
+            ;;
+        (-decompose-dict=*)
+            optValue="${1#*=}"
+            case "$optValue" in
+            ('' | none | false) ;;  ## Ignore
+            (*)
+                appArgs="$appArgs -decomposeParDict $optValue"
+                nProcs="$(getNumberOfProcessors "$optValue")"
                 ;;
-            -o | -overwrite)
-                logMode=overwrite
-                ;;
-            -s | -suffix)
-                logFile=".$2"
-                shift
-                ;;
-            -n | -np)
-                nProcs="$2"
-                shift
-                ;;
-            -decomposeParDict)
-                appArgs="$appArgs $1 $2"
-                nProcs=$(getNumberOfProcessors "$2")
-                shift
-                ;;
-            '')
-                ;;
-            *)
-                appRun="$1"
+            esac
+            ;;
+        (-decomposeParDict)
+            optValue="$2"
+            shift
+            case "$optValue" in
+            ('' | none | false) ;;  ## Ignore
+            (*)
+                appArgs="$appArgs -decomposeParDict $optValue"
+                nProcs="$(getNumberOfProcessors "$optValue")"
                 ;;
+            esac
+            ;;
+         (*)
+            appRun="$1"
+            ;;
         esac
         shift
     done

src/OpenFOAM/containers/PtrLists/PtrList/PtrList.H

@@ -173,6 +173,10 @@ public:
         //- (discard old element at that location).
         //- Return reference to the new list element.
         template<class... Args>
+        inline T& emplace_set(const label i, Args&&... args);
+
+        //- Same as emplace_set()
+        template<class... Args>
         inline T& emplace(const label i, Args&&... args);
 
         //- Set element to given pointer and return old element (can be null)

src/OpenFOAM/containers/PtrLists/PtrList/PtrListI.H

@@ -172,11 +172,17 @@ inline void Foam::PtrList<T>::push_back(PtrList<T>&& other)
 
 template<class T>
 template<class... Args>
-inline T& Foam::PtrList<T>::emplace
-(
-    const label i,
-    Args&&... args
-)
+inline T& Foam::PtrList<T>::emplace_set(const label i, Args&&... args)
+{
+    T* ptr = new T(std::forward<Args>(args)...);
+    (void)this->set(i, ptr);
+    return *ptr;
+}
+
+
+template<class T>
+template<class... Args>
+inline T& Foam::PtrList<T>::emplace(const label i, Args&&... args)
 {
     T* ptr = new T(std::forward<Args>(args)...);
     (void)this->set(i, ptr);

src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H

@@ -470,8 +470,7 @@ public:
 
     // Exchange
 
-        //- Helper: exchange sizes of sendBufs for specified
-        //- set of send/receive processes.
+        //- Helper: exchange sizes of sendBufs for specified send/recv ranks
         template<class Container>
         static void exchangeSizes
         (
@@ -483,6 +482,17 @@ public:
             const label comm = UPstream::worldComm
         );
 
+        //- Helper: exchange sizes of sendBufs for specified neighbour ranks
+        template<class Container>
+        static void exchangeSizes
+        (
+            const labelUList& neighProcs,
+            const Container& sendBufs,
+            labelList& sizes,
+            const label tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm
+        );
+
         //- Helper: exchange sizes of sendBufs.
         //- The sendBufs is the data per processor (in the communicator).
         //  Returns sizes of sendBufs on the sending processor.

src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C

@@ -28,6 +28,30 @@ License
 
 #include "PstreamBuffers.H"
 #include "bitSet.H"
+#include "debug.H"
+#include "registerSwitch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::PstreamBuffers::algorithm
+(
+    // Not really the most creative name...
+    Foam::debug::optimisationSwitch("pbufs.algorithm", -1)
+);
+registerOptSwitch
+(
+    "pbufs.algorithm",
+    int,
+    Foam::PstreamBuffers::algorithm
+);
+
+
+// Simple enumerations
+// -------------------
+static constexpr int algorithm_PEX_allToAll = -1;  // OpenFOAM 2212 and earlier
+//static constexpr int algorithm_PEX_hybrid = 0;   // New default?
+static constexpr int algorithm_full_NBX = 1;       // Experimental
+
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -43,9 +67,58 @@ void Foam::PstreamBuffers::finalExchange
 
     if (commsType_ == UPstream::commsTypes::nonBlocking)
     {
-        // Dense storage uses all-to-all
-        Pstream::exchangeSizes(sendBuffers_, recvSizes, comm_);
+        if
+        (
+            wait
+         && (algorithm >= algorithm_full_NBX)
+         && (UPstream::maxCommsSize <= 0)
+        )
+        {
+            // NBX algorithm (nonblocking exchange)
+            // - when requested and waiting, no data chunking etc
 
+            PstreamDetail::exchangeConsensus<DynamicList<char>, char>
+            (
+                sendBuffers_,
+                recvBuffers_,
+                recvSizes,
+                (tag_ + 271828),  // some unique tag?
+                comm_,
+                wait
+            );
+
+            return;
+        }
+
+
+        // PEX algorithm with two different flavours of exchanging sizes
+
+        // Assemble the send sizes (cf. Pstream::exchangeSizes)
+        labelList sendSizes(nProcs_);
+        forAll(sendBuffers_, proci)
+        {
+            sendSizes[proci] = sendBuffers_[proci].size();
+        }
+        recvSizes.resize_nocopy(nProcs_);
+
+        if (algorithm == algorithm_PEX_allToAll)
+        {
+            // PEX stage 1: exchange sizes (all-to-all)
+            UPstream::allToAll(sendSizes, recvSizes, comm_);
+        }
+        else
+        {
+            // PEX stage 1: exchange sizes (non-blocking consensus)
+            UPstream::allToAllConsensus
+            (
+                sendSizes,
+                recvSizes,
+                (tag_ + 314159),  // some unique tag?
+                comm_
+            );
+        }
+
+        // PEX stage 2: point-to-point data exchange
         Pstream::exchange<DynamicList<char>, char>
         (
             sendBuffers_,
@@ -73,6 +146,30 @@ void Foam::PstreamBuffers::finalExchange
 
     if (commsType_ == UPstream::commsTypes::nonBlocking)
     {
+        // Preparation. Temporarily abuse recvSizes as logic to clear
+        // send buffers that are not in the neighbourhood connection
+        {
+            recvSizes.resize_nocopy(nProcs_);
+            recvSizes = 0;
+
+            // Preserve self-send, even if not described by neighbourhood
+            recvSizes[UPstream::myProcNo(comm_)] = 1;
+
+            for (const label proci : sendProcs)
+            {
+                recvSizes[proci] = 1;  // Connected
+            }
+
+            for (label proci=0; proci < nProcs_; ++proci)
+            {
+                if (!recvSizes[proci])  // Not connected
+                {
+                    sendBuffers_[proci].clear();
+                }
+            }
+        }
+
+        // PEX stage 1: exchange sizes (limited neighbourhood)
         Pstream::exchangeSizes
         (
             sendProcs,
@@ -83,6 +180,7 @@ void Foam::PstreamBuffers::finalExchange
             comm_
         );
 
+        // PEX stage 2: point-to-point data exchange
         Pstream::exchange<DynamicList<char>, char>
         (
             sendBuffers_,
@@ -96,11 +194,10 @@ void Foam::PstreamBuffers::finalExchange
 }
 
 
-void Foam::PstreamBuffers::finalExchangeGatherScatter
+void Foam::PstreamBuffers::finalGatherScatter
 (
     const bool isGather,
     const bool wait,
-    const bool needSizes,
     labelList& recvSizes
 )
 {
@@ -133,6 +230,12 @@ void Foam::PstreamBuffers::finalExchangeGatherScatter
 
     if (commsType_ == UPstream::commsTypes::nonBlocking)
     {
+        // Use PEX algorithm
+        // - for a non-sparse gather/scatter, it is presumed that
+        // MPI_Gather/MPI_Scatter will be the most efficient way to
+        // communicate the sizes.
+
+        // PEX stage 1: exchange sizes (gather or scatter)
         if (isGather)
         {
             // gather mode (all-to-one): master [0] <- everyone
@@ -166,7 +269,7 @@ void Foam::PstreamBuffers::finalExchangeGatherScatter
             recvSizes[0] = myRecv;
         }
 
-
+        // PEX stage 2: point-to-point data exchange
         Pstream::exchange<DynamicList<char>, char>
         (
             sendBuffers_,
@@ -378,7 +481,7 @@ Foam::label Foam::PstreamBuffers::recvDataCount(const label proci) const
 
 Foam::labelList Foam::PstreamBuffers::recvDataCounts() const
 {
-    labelList counts(recvPositions_.size(), Zero);
+    labelList counts(nProcs_, Zero);
 
     if (finishedSendsCalled_)
     {
@@ -404,17 +507,20 @@ Foam::labelList Foam::PstreamBuffers::recvDataCounts() const
 }
 
 
-Foam::label Foam::PstreamBuffers::maxNonLocalRecvCount(const label proci) const
+Foam::label Foam::PstreamBuffers::maxNonLocalRecvCount
+(
+    const label excludeProci
+) const
 {
     label maxLen = 0;
 
     if (finishedSendsCalled_)
     {
-        forAll(recvBuffers_, idx)
+        forAll(recvBuffers_, proci)
         {
-            const label len(recvBuffers_[idx].size() - recvPositions_[idx]);
-            if (idx != proci)
+            if (excludeProci != proci)
             {
+                label len(recvBuffers_[proci].size() - recvPositions_[proci]);
                 maxLen = max(maxLen, len);
             }
         }
@@ -494,6 +600,9 @@ void Foam::PstreamBuffers::finishedSends
     const bool wait
 )
 {
+    // Resize for copying back
+    recvSizes.resize_nocopy(sendBuffers_.size());
+
     finalExchange(wait, recvSizes);
 
     if (commsType_ != UPstream::commsTypes::nonBlocking)
@@ -510,39 +619,25 @@ void Foam::PstreamBuffers::finishedSends
 }
 
 
-void Foam::PstreamBuffers::finishedSends
+void Foam::PstreamBuffers::finishedNeighbourSends
 (
-    const labelUList& sendProcs,
-    const labelUList& recvProcs,
-    const bool wait
-)
-{
-    labelList recvSizes;
-    finalExchange(sendProcs, recvProcs, wait, recvSizes);
-}
-
-
-void Foam::PstreamBuffers::finishedSends
-(
-    const labelUList& sendProcs,
-    const labelUList& recvProcs,
+    const labelUList& neighProcs,
     labelList& recvSizes,
     const bool wait
 )
 {
-    finalExchange(sendProcs, recvProcs, wait, recvSizes);
+    finalExchange(neighProcs, neighProcs, wait, recvSizes);
+}
 
-    if (commsType_ != UPstream::commsTypes::nonBlocking)
-    {
-        FatalErrorInFunction
-            << "Obtaining sizes not supported in "
-            << UPstream::commsTypeNames[commsType_] << endl
-            << " since transfers already in progress. Use non-blocking instead."
-            << exit(FatalError);
 
-        // Note: maybe possible only if using different tag from write started
-        // by ~UOPstream. Needs some work.
-    }
+void Foam::PstreamBuffers::finishedNeighbourSends
+(
+    const labelUList& neighProcs,
+    const bool wait
+)
+{
+    labelList recvSizes;
+    finalExchange(neighProcs, neighProcs, wait, recvSizes);
 }
 
 
@@ -562,10 +657,8 @@ bool Foam::PstreamBuffers::finishedSends
     }
 
     // Update send connections
-    // - reasonable to assume there are no self-sends on UPstream::myProcNo
     forAll(sendBuffers_, proci)
     {
-        // ie, sendDataCount(proci) != 0
         if (sendConnections.set(proci, !sendBuffers_[proci].empty()))
         {
             // The state changed
@@ -577,22 +670,20 @@ bool Foam::PstreamBuffers::finishedSends
 
     if (changed)
     {
-        // Create send/recv topology
+        // Update send/recv topology
+        labelList recvSizes;
+        finishedSends(recvSizes, wait);  // eg, using all-to-all
 
         // The send ranks
         sendProcs.clear();
         forAll(sendBuffers_, proci)
         {
-            // ie, sendDataCount(proci) != 0
             if (!sendBuffers_[proci].empty())
             {
                 sendProcs.push_back(proci);
             }
         }
 
-        labelList recvSizes;
-        finishedSends(recvSizes, wait);  // All-to-all
-
         // The recv ranks
         recvProcs.clear();
         forAll(recvSizes, proci)
@@ -606,46 +697,25 @@ bool Foam::PstreamBuffers::finishedSends
     else
     {
         // Use existing send/recv ranks
-
-        finishedSends(sendProcs, recvProcs, wait);
+        labelList recvSizes;
+        finalExchange(sendProcs, recvProcs, wait, recvSizes);
     }
 
     return changed;
 }
 
 
-void Foam::PstreamBuffers::finishedNeighbourSends
-(
-    const labelUList& neighProcs,
-    labelList& recvSizes,
-    const bool wait
-)
-{
-    finishedSends(neighProcs, neighProcs, recvSizes, wait);
-}
-
-
-void Foam::PstreamBuffers::finishedNeighbourSends
-(
-    const labelUList& neighProcs,
-    const bool wait
-)
-{
-    finishedSends(neighProcs, neighProcs, wait);
-}
-
-
 void Foam::PstreamBuffers::finishedGathers(const bool wait)
 {
     labelList recvSizes;
-    finalExchangeGatherScatter(true, wait, false, recvSizes);
+    finalGatherScatter(true, wait, recvSizes);
 }
 
 
 void Foam::PstreamBuffers::finishedScatters(const bool wait)
 {
     labelList recvSizes;
-    finalExchangeGatherScatter(false, wait, false, recvSizes);
+    finalGatherScatter(false, wait, recvSizes);
 }
 
 
@@ -655,7 +725,7 @@ void Foam::PstreamBuffers::finishedGathers
     const bool wait
 )
 {
-    finalExchangeGatherScatter(true, wait, true, recvSizes);
+    finalGatherScatter(true, wait, recvSizes);
 
     if (commsType_ != UPstream::commsTypes::nonBlocking)
     {
@@ -677,7 +747,7 @@ void Foam::PstreamBuffers::finishedScatters
     const bool wait
 )
 {
-    finalExchangeGatherScatter(false, wait, true, recvSizes);
+    finalGatherScatter(false, wait, recvSizes);
 
     if (commsType_ != UPstream::commsTypes::nonBlocking)
     {

src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H

@@ -164,7 +164,7 @@ class PstreamBuffers
         );
 
         //- Mark sends as done.
-        //  Only exchange sizes using the sendProcs/recvProcs subset
+        //  Only exchange sizes using the neighbour ranks
         //  (nonBlocking comms).
         void finalExchange
         (
@@ -175,11 +175,10 @@ class PstreamBuffers
         );
 
         //- For all-to-one or one-to-all
-        void finalExchangeGatherScatter
+        void finalGatherScatter
         (
             const bool isGather,
             const bool wait,
-            const bool needSizes,   // If recvSizes needed or scratch
             labelList& recvSizes
         );
 
@@ -202,6 +201,12 @@ class PstreamBuffers
 
 public:
 
+    // Static Data
+
+        //- Preferred exchange algorithm (may change or be removed in future)
+        static int algorithm;
+
+
     // Constructors
 
         //- Construct given communication type (default: nonBlocking), message
@@ -340,7 +345,7 @@ public:
 
         //- Maximum receive size, excluding the specified processor rank
         //- Must call finishedSends() or other finished.. method first!
-        label maxNonLocalRecvCount(const label proci) const;
+        label maxNonLocalRecvCount(const label excludeProci) const;
 
         //- Number of unconsumed receive bytes for the specified processor.
         //- Must call finishedSends() or other finished.. method first!
@@ -397,60 +402,6 @@ public:
 
     // Functions with restricted neighbours
 
-        //- Mark sends as done using subset of send/recv ranks
-        //- to exchange data on.
-        //
-        //  Non-blocking mode: populates receive buffers.
-        //  \param sendProcs ranks used for sends
-        //  \param recvProcs ranks used for recvs
-        //  \param wait wait for requests to complete (in nonBlocking mode)
-        //
-        //  \warning currently only valid for nonBlocking comms.
-        void finishedSends
-        (
-            const labelUList& sendProcs,
-            const labelUList& recvProcs,
-            const bool wait = true
-        );
-
-        //- Mark sends as done using subset of send/recv ranks
-        //- to exchange data on. Recovers the sizes (bytes) received.
-        //
-        //  Non-blocking mode: populates receive buffers.
-        //
-        //  \param sendProcs ranks used for sends
-        //  \param recvProcs ranks used for recvs
-        //  \param[out] recvSizes the sizes (bytes) received
-        //  \param wait wait for requests to complete (in nonBlocking mode)
-        //
-        //  \warning currently only valid for nonBlocking comms.
-        void finishedSends
-        (
-            const labelUList& sendProcs,
-            const labelUList& recvProcs,
-            labelList& recvSizes,
-            const bool wait = true
-        );
-
-        //- A caching version that uses a limited send/recv connectivity.
-        //
-        //  Non-blocking mode: populates receive buffers.
-        //  \param sendConnections on/off for sending ranks
-        //  \param sendProcs ranks used for sends
-        //  \param recvProcs ranks used for recvs
-        //  \param wait wait for requests to complete (in nonBlocking mode)
-        //
-        //  \return True if the send/recv connectivity changed
-        //
-        //  \warning currently only valid for nonBlocking comms.
-        bool finishedSends
-        (
-            bitSet& sendConnections,
-            DynamicList<label>& sendProcs,
-            DynamicList<label>& recvProcs,
-            const bool wait = true
-        );
-
         //- Mark sends as done using subset of send/recv ranks
         //- and recover the sizes (bytes) received.
         //
@@ -484,6 +435,25 @@ public:
             const bool wait = true
         );
 
+        //- A caching version that uses a limited send/recv connectivity.
+        //
+        //  Non-blocking mode: populates receive buffers.
+        //  \param sendConnections on/off for sending ranks
+        //  \param sendProcs ranks used for sends
+        //  \param recvProcs ranks used for recvs
+        //  \param wait wait for requests to complete (in nonBlocking mode)
+        //
+        //  \return True if the send/recv connectivity changed
+        //
+        //  \warning currently only valid for nonBlocking comms.
+        bool finishedSends
+        (
+            bitSet& sendConnections,
+            DynamicList<label>& sendProcs,
+            DynamicList<label>& recvProcs,
+            const bool wait = true
+        );
+
 
     // Gather/scatter modes
 

src/OpenFOAM/db/IOstreams/Pstreams/PstreamExchange.C

@@ -746,58 +746,89 @@ void Foam::Pstream::exchangeSizes
     const label comm
 )
 {
-    //const label myProci = UPstream::myProcNo(comm);
+    const label myProci = UPstream::myProcNo(comm);
+    const label numProcs = UPstream::nProcs(comm);
 
-    if (sendBufs.size() != UPstream::nProcs(comm))
+    if (sendBufs.size() != numProcs)
     {
         FatalErrorInFunction
             << "Size of container " << sendBufs.size()
-            << " does not equal the number of processors "
-            << UPstream::nProcs(comm)
+            << " does not equal the number of processors " << numProcs
             << Foam::abort(FatalError);
     }
 
-    labelList sendSizes(sendProcs.size());
-    forAll(sendProcs, i)
+    labelList sendSizes(numProcs);
+    for (label proci = 0; proci < numProcs; ++proci)
     {
-        sendSizes[i] = sendBufs[sendProcs[i]].size();
+        sendSizes[proci] = sendBufs[proci].size();
     }
 
-    recvSizes.resize_nocopy(sendBufs.size());
+    recvSizes.resize_nocopy(numProcs);
     recvSizes = 0;  // Ensure non-received entries are properly zeroed
 
+    // Preserve self-send, even if not described by neighbourhood
+    recvSizes[myProci] = sendSizes[myProci];
+
     const label startOfRequests = UPstream::nRequests();
 
     for (const label proci : recvProcs)
     {
-        UIPstream::read
-        (
-            UPstream::commsTypes::nonBlocking,
-            proci,
-            reinterpret_cast<char*>(&recvSizes[proci]),
-            sizeof(label),
-            tag,
-            comm
-        );
+        if (proci != myProci)
+        {
+            UIPstream::read
+            (
+                UPstream::commsTypes::nonBlocking,
+                proci,
+                reinterpret_cast<char*>(&recvSizes[proci]),
+                sizeof(label),
+                tag,
+                comm
+            );
+        }
     }
 
-    forAll(sendProcs, i)
+    for (const label proci : sendProcs)
     {
-        UOPstream::write
-        (
-            UPstream::commsTypes::nonBlocking,
-            sendProcs[i],
-            reinterpret_cast<char*>(&sendSizes[i]),
-            sizeof(label),
-            tag,
-            comm
-        );
+        if (proci != myProci)
+        {
+            UOPstream::write
+            (
+                UPstream::commsTypes::nonBlocking,
+                proci,
+                reinterpret_cast<char*>(&sendSizes[proci]),
+                sizeof(label),
+                tag,
+                comm
+            );
+        }
     }
 
     UPstream::waitRequests(startOfRequests);
 }
 
 
+template<class Container>
+void Foam::Pstream::exchangeSizes
+(
+    const labelUList& neighProcs,
+    const Container& sendBufs,
+    labelList& recvSizes,
+    const label tag,
+    const label comm
+)
+{
+    Pstream::exchangeSizes<Container>
+    (
+        neighProcs,  // send
+        neighProcs,  // recv
+        sendBufs,
+        recvSizes,
+        tag,
+        comm
+    );
+}
+
+
 // Sparse sending
 template<class Container>
 void Foam::Pstream::exchangeSizes

src/OpenFOAM/global/argList/argList.C

@@ -1720,6 +1720,8 @@ void Foam::argList::parse
                 << "    nProcsSimpleSum    : " << Pstream::nProcsSimpleSum << nl
                 << "    nonBlockingExchange: "
                 << Pstream::nProcsNonblockingExchange << nl
+                << "    exchange algorithm : "
+                << PstreamBuffers::algorithm << nl
                 << "    commsType          : "
                 << Pstream::commsTypeNames[Pstream::defaultCommsType] << nl
                 << "    polling iterations : " << Pstream::nPollProcInterfaces

src/OpenFOAM/meshes/boundBox/boundBox.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -163,6 +163,9 @@ public:
         //- Construct a bounding box containing a single initial point
         inline explicit boundBox(const point& p);
 
+        //- Construct a 0/1 unit bounding box
+        inline explicit boundBox(const Foam::zero_one);
+
         //- Construct from bound box min/max points
         inline boundBox(const point& min, const point& max);
 
@@ -291,6 +294,9 @@ public:
             //- Reset to an inverted box
             inline void reset();
 
+            //- Reset to a 0/1 unit bounding box
+            inline void reset(const Foam::zero_one);
+
             //- Reset min/max to be identical to the specified point
             inline void reset(const point& pt);
 

src/OpenFOAM/meshes/boundBox/boundBoxI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -112,6 +112,13 @@ inline Foam::boundBox::boundBox()
 {}
 
 
+inline Foam::boundBox::boundBox(const Foam::zero_one)
+:
+    min_(point::zero),
+    max_(point::one)
+{}
+
+
 inline Foam::boundBox::boundBox(const point& p)
 :
     min_(p),
@@ -299,6 +306,13 @@ inline void Foam::boundBox::reset()
 }
 
 
+inline void Foam::boundBox::reset(const Foam::zero_one)
+{
+    min_ = point::zero;
+    max_ = point::one;
+}
+
+
 inline void Foam::boundBox::reset(const point& pt)
 {
     min_ = pt;

src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ijkAddressing_H
-#define ijkAddressing_H
+#ifndef Foam_ijkAddressing_H
+#define Foam_ijkAddressing_H
 
 #include "labelVector.H"
 #include "vector.H"

src/OpenFOAM/meshes/ijkMesh/ijkMesh.H

@@ -38,8 +38,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef ijkMesh_H
-#define ijkMesh_H
+#ifndef Foam_ijkMesh_H
+#define Foam_ijkMesh_H
 
 #include "ijkAddressing.H"
 

src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -201,6 +201,9 @@ public:
 
     // Constructors
 
+        //- Construct a 0/1 unit bounding box
+        inline explicit treeBoundBox(const Foam::zero_one);
+
         //- Construct a bounding box containing a single initial point
         inline explicit treeBoundBox(const point& p);
 

src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,6 +31,12 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+inline Foam::treeBoundBox::treeBoundBox(const Foam::zero_one)
+:
+    boundBox(zero_one{})
+{}
+
+
 inline Foam::treeBoundBox::treeBoundBox(const point& p)
 :
     boundBox(p)

src/OpenFOAM/parallel/globalIndex/globalIndex.C

@@ -214,7 +214,7 @@ void Foam::globalIndex::reset
             // Non-parallel branch: use localSize on-proc, zero elsewhere
 
             localLens.resize(len, Zero);
-            localLens[Pstream::myProcNo(comm)] = localSize;
+            localLens[UPstream::myProcNo(comm)] = localSize;
         }
 
         reset(localLens, true);  // checkOverflow = true

src/OpenFOAM/parallel/globalIndex/globalIndex.H

@@ -519,7 +519,9 @@ public:
         {
             List<Type> allData;
             gather(offsets, comm, procIDs, fld, allData, tag, ct);
-            if (Pstream::myProcNo(comm) == procIDs[0])
+
+            const int masterProci = procIDs.size() ? procIDs[0] : 0;
+            if (UPstream::myProcNo(comm) == masterProci)
             {
                 fld.transfer(allData);
             }

src/OpenFOAM/parallel/globalIndex/globalIndexI.H

@@ -137,14 +137,14 @@ inline Foam::globalIndex::globalIndex
 
 inline bool Foam::globalIndex::empty() const
 {
-    return offsets_.empty() || offsets_.last() == 0;
+    return offsets_.empty() || offsets_.back() == 0;
 }
 
 
 inline Foam::label Foam::globalIndex::totalSize() const
 {
     const label len = (offsets_.size() - 1);
-    return (len < 1) ? static_cast<label>(0) : offsets_[len];
+    return (len < 1) ? 0 : offsets_[len];
 }
 
 
@@ -163,7 +163,7 @@ inline Foam::labelList Foam::globalIndex::sizes() const
 inline Foam::label Foam::globalIndex::nProcs() const noexcept
 {
     const label len = (offsets_.size() - 1);
-    return (len < 1) ? static_cast<label>(0) : len;
+    return (len < 1) ? 0 : len;
 }
 
 
@@ -219,7 +219,7 @@ inline Foam::label Foam::globalIndex::localStart(const label proci) const
 
 inline Foam::label Foam::globalIndex::localStart() const
 {
-    return localStart(Pstream::myProcNo());
+    return localStart(UPstream::myProcNo());
 }
 
 
@@ -231,7 +231,7 @@ inline Foam::label Foam::globalIndex::localSize(const label proci) const
 
 inline Foam::label Foam::globalIndex::localSize() const
 {
-    return localSize(Pstream::myProcNo());
+    return localSize(UPstream::myProcNo());
 }
 
 
@@ -244,7 +244,7 @@ inline Foam::label Foam::globalIndex::maxSize() const
 
 inline Foam::label Foam::globalIndex::maxNonLocalSize() const
 {
-    return maxNonLocalSize(Pstream::myProcNo());
+    return maxNonLocalSize(UPstream::myProcNo());
 }
 
 
@@ -256,7 +256,7 @@ inline Foam::labelRange Foam::globalIndex::range(const label proci) const
 
 inline Foam::labelRange Foam::globalIndex::range() const
 {
-    return range(Pstream::myProcNo());
+    return range(UPstream::myProcNo());
 }
 
 
@@ -268,7 +268,7 @@ inline bool Foam::globalIndex::isLocal(const label proci, const label i) const
 
 inline bool Foam::globalIndex::isLocal(const label i) const
 {
-    return isLocal(Pstream::myProcNo(), i);
+    return isLocal(UPstream::myProcNo(), i);
 }
 
 
@@ -284,7 +284,7 @@ inline Foam::label Foam::globalIndex::toGlobal
 
 inline Foam::label Foam::globalIndex::toGlobal(const label i) const
 {
-    return toGlobal(Pstream::myProcNo(), i);
+    return toGlobal(UPstream::myProcNo(), i);
 }
 
 
@@ -306,7 +306,7 @@ inline Foam::labelList Foam::globalIndex::toGlobal
     const labelUList& labels
 ) const
 {
-    return toGlobal(Pstream::myProcNo(), labels);
+    return toGlobal(UPstream::myProcNo(), labels);
 }
 
 
@@ -327,7 +327,7 @@ inline void Foam::globalIndex::inplaceToGlobal
 
 inline void Foam::globalIndex::inplaceToGlobal(labelUList& labels) const
 {
-    inplaceToGlobal(Pstream::myProcNo(), labels);
+    inplaceToGlobal(UPstream::myProcNo(), labels);
 }
 
 
@@ -349,7 +349,7 @@ Foam::globalIndex::toLocal(const label proci, const label i) const
 
 inline Foam::label Foam::globalIndex::toLocal(const label i) const
 {
-    return toLocal(Pstream::myProcNo(), i);
+    return toLocal(UPstream::myProcNo(), i);
 }
 
 
@@ -363,7 +363,7 @@ inline Foam::label Foam::globalIndex::whichProcID(const label i) const
             << abort(FatalError);
     }
 
-    const label proci(Pstream::myProcNo());
+    const int proci = UPstream::myProcNo();
 
     return isLocal(proci, i) ? proci : findLower(offsets_, i+1);
 }

src/OpenFOAM/parallel/globalIndex/globalIndexTemplates.C

@@ -91,7 +91,9 @@ void Foam::globalIndex::gatherValues
 
     const label startOfRequests = UPstream::nRequests();
 
-    if (UPstream::myProcNo(comm) == procIDs[0])
+    const int masterProci = procIDs.size() ? procIDs[0] : 0;
+
+    if (UPstream::myProcNo(comm) == masterProci)
     {
         allValues.resize_nocopy(procIDs.size());
         allValues[0] = localValue;
@@ -126,7 +128,7 @@ void Foam::globalIndex::gatherValues
             UOPstream::write
             (
                 commsType,
-                procIDs[0],
+                masterProci,
                 reinterpret_cast<const char*>(&localValue),
                 sizeof(Type),
                 tag,
@@ -135,7 +137,7 @@ void Foam::globalIndex::gatherValues
         }
         else
         {
-            OPstream toMaster(commsType, procIDs[0], 0, tag, comm);
+            OPstream toMaster(commsType, masterProci, 0, tag, comm);
             toMaster << localValue;
         }
     }
@@ -176,9 +178,11 @@ void Foam::globalIndex::gather
 
     const label startOfRequests = UPstream::nRequests();
 
-    if (Pstream::myProcNo(comm) == procIDs[0])
+    const int masterProci = procIDs.size() ? procIDs[0] : 0;
+
+    if (UPstream::myProcNo(comm) == masterProci)
     {
-        allFld.resize_nocopy(off.last());  // == totalSize()
+        allFld.resize_nocopy(off.back());  // == totalSize()
 
         // Assign my local data - respect offset information
         // so that we can request 0 entries to be copied.
@@ -226,7 +230,7 @@ void Foam::globalIndex::gather
             UOPstream::write
             (
                 commsType,
-                procIDs[0],
+                masterProci,
                 fld.cdata_bytes(),
                 fld.size_bytes(),
                 tag,
@@ -235,7 +239,7 @@ void Foam::globalIndex::gather
         }
         else
         {
-            OPstream toMaster(commsType, procIDs[0], 0, tag, comm);
+            OPstream toMaster(commsType, masterProci, 0, tag, comm);
             toMaster << fld;
         }
     }
@@ -294,9 +298,11 @@ void Foam::globalIndex::gather
 
     const label startOfRequests = UPstream::nRequests();
 
-    if (Pstream::myProcNo(comm) == procIDs[0])
+    const int masterProci = procIDs.size() ? procIDs[0] : 0;
+
+    if (UPstream::myProcNo(comm) == masterProci)
     {
-        allFld.resize_nocopy(off.last());  // == totalSize()
+        allFld.resize_nocopy(off.back());  // == totalSize()
 
         // Assign my local data - respect offset information
         // so that we can request 0 entries to be copied
@@ -331,7 +337,7 @@ void Foam::globalIndex::gather
         }
         else
         {
-            OPstream toMaster(commsType, procIDs[0], 0, tag, comm);
+            OPstream toMaster(commsType, masterProci, 0, tag, comm);
             toMaster << fld;
         }
     }
@@ -877,7 +883,9 @@ void Foam::globalIndex::scatter
 
     const label startOfRequests = UPstream::nRequests();
 
-    if (Pstream::myProcNo(comm) == procIDs[0])
+    const int masterProci = procIDs.size() ? procIDs[0] : 0;
+
+    if (UPstream::myProcNo(comm) == masterProci)
     {
         for (label i = 1; i < procIDs.size(); ++i)
         {
@@ -929,7 +937,7 @@ void Foam::globalIndex::scatter
             UIPstream::read
             (
                 commsType,
-                procIDs[0],
+                masterProci,
                 fld.data_bytes(),
                 fld.size_bytes(),
                 tag,
@@ -938,7 +946,7 @@ void Foam::globalIndex::scatter
         }
         else
         {
-            IPstream fromMaster(commsType, procIDs[0], 0, tag, comm);
+            IPstream fromMaster(commsType, masterProci, 0, tag, comm);
             fromMaster >> fld;
         }
     }

src/fileFormats/vtk/file/foamVtkFileWriter.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -76,11 +76,11 @@ Foam::Ostream& Foam::vtk::fileWriter::reportBadState
 Foam::Ostream& Foam::vtk::fileWriter::reportBadState
 (
     Ostream& os,
-    outputState expected,
+    outputState expected1,
     outputState expected2
 ) const
 {
-    reportBadState(os, expected)
+    reportBadState(os, expected1)
         << " or (" << stateNames[expected2] << ')';
     return os;
 }
@@ -554,10 +554,18 @@ bool Foam::vtk::fileWriter::writeProcIDs(const label nValues)
     {
         ++nCellData_;
     }
+    else if (isState(outputState::POINT_DATA))
+    {
+        ++nPointData_;
+    }
     else
     {
-        reportBadState(FatalErrorInFunction, outputState::CELL_DATA)
-            << " for procID field" << nl << endl
+        reportBadState
+        (
+            FatalErrorInFunction,
+            outputState::CELL_DATA,
+            outputState::POINT_DATA
+        )   << " for procID field" << nl << endl
             << exit(FatalError);
 
         return false;

src/fileFormats/vtk/file/foamVtkFileWriter.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -201,7 +201,8 @@ protected:
             const UList<Type>& field
         );
 
-        //- Write nValues of processor ids as CellData (no-op in serial)
+        //- Write nValues of processor ids as CellData or PointData
+        //- (no-op in serial)
         bool writeProcIDs(const label nValues);
 
 
@@ -305,10 +306,16 @@ public:
         //  \return True if the state changed
         virtual bool beginPointData(label nFields = 0) = 0;
 
-        //- Return the number of CellData written for the Piece thus far.
+        //- True if output state corresponds to CELL_DATA
+        inline bool isCellData() const noexcept;
+
+        //- True if output state corresponds to POINT_DATA
+        inline bool isPointData() const noexcept;
+
+        //- The number of CellData written for the Piece thus far.
         inline label nCellData() const noexcept;
 
-        //- Return the number of PointData written for the Piece thus far.
+        //- The number of PointData written for the Piece thus far.
         inline label nPointData() const noexcept;
 
 

src/fileFormats/vtk/file/foamVtkFileWriterI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -107,4 +107,16 @@ inline Foam::label Foam::vtk::fileWriter::nPointData() const noexcept
 }
 
 
+inline bool Foam::vtk::fileWriter::isCellData() const noexcept
+{
+    return (outputState::CELL_DATA == state_);
+}
+
+
+inline bool Foam::vtk::fileWriter::isPointData() const noexcept
+{
+    return (outputState::POINT_DATA == state_);
+}
+
+
 // ************************************************************************* //

src/fileFormats/vtk/write/foamVtkLineWriter.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -131,6 +131,10 @@ void Foam::vtk::lineWriter::piece
 
 bool Foam::vtk::lineWriter::writeProcIDs()
 {
+    if (this->isPointData())
+    {
+        return vtk::fileWriter::writeProcIDs(nLocalPoints_);
+    }
     return vtk::fileWriter::writeProcIDs(nLocalLines_);
 }
 

src/fileFormats/vtk/write/foamVtkLineWriter.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -155,8 +155,8 @@ public:
         void piece(const pointField& points, const edgeList& edges);
 
 
-        //- Write processor ids for each line as CellData
-        //- (no-op in serial)
+        //- Write processor ids for each line as CellData or for each point
+        //- as PointData, depending on isPointData() state. No-op in serial.
         bool writeProcIDs();
 
         //- Write a uniform field of Cell (Line) or Point values

src/fileFormats/vtk/write/foamVtkSurfaceWriter.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -131,6 +131,10 @@ void Foam::vtk::surfaceWriter::piece
 
 bool Foam::vtk::surfaceWriter::writeProcIDs()
 {
+    if (this->isPointData())
+    {
+        return vtk::fileWriter::writeProcIDs(nLocalPoints_);
+    }
     return vtk::fileWriter::writeProcIDs(nLocalPolys_);
 }
 

src/fileFormats/vtk/write/foamVtkSurfaceWriter.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -155,8 +155,8 @@ public:
         void piece(const pointField& points, const faceList& faces);
 
 
-        //- Write processor ids for each poly as CellData
-        //- (no-op in serial)
+        //- Write processor ids for each poly as CellData or for each point
+        //- as PointData, depending on isPointData() state. No-op in serial.
         bool writeProcIDs();
 
         //- Write a uniform field of Cell (Poly) or Point values

src/finiteVolume/fvMesh/zoneDistribute/zoneDistribute.C

@@ -114,13 +114,15 @@ void Foam::zoneDistribute::setUpCommforZone
 
         pBufs_.clear();
 
-        for (const int proci : UPstream::allProcs())
+        for (const int proci : pBufs_.allProcs())
         {
-            if (proci != UPstream::myProcNo() && !needed[proci].empty())
+            const auto& indices = needed[proci];
+
+            if (proci != UPstream::myProcNo() && !indices.empty())
             {
                 // Serialize as List
                 UOPstream toProc(proci, pBufs_);
-                toProc << needed[proci].sortedToc();
+                toProc << indices.sortedToc();
             }
         }
 

src/finiteVolume/fvMesh/zoneDistribute/zoneDistributeI.H

@@ -175,14 +175,16 @@ Foam::Map<Type> Foam::zoneDistribute::getDatafromOtherProc
 
         pBufs_.clear();
 
-        for (const int proci : UPstream::allProcs())
+        for (const int proci : pBufs_.allProcs())
         {
-            if (proci != UPstream::myProcNo() && !send_[proci].empty())
+            const auto& indices = send_[proci];
+
+            if (proci != UPstream::myProcNo() && !indices.empty())
             {
                 // Serialize as Map
-                Map<Type> sendValues(2*send_[proci].size());
+                Map<Type> sendValues(2*indices.size());
 
-                for (const label sendIdx : send_[proci])
+                for (const label sendIdx : indices)
                 {
                     sendValues.insert
                     (

src/mesh/blockMesh/PDRblockMesh/PDRblock.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -114,6 +114,22 @@ const Foam::PDRblock& Foam::PDRblock::null()
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::PDRblock::addDefaultPatches()
+{
+    // Default boundaries with patchi == shapeFacei
+    patches_.resize(6);
+    for (label patchi=0; patchi < 6; ++patchi)
+    {
+        boundaryEntry& bentry = patches_.emplace_set(patchi);
+
+        bentry.name_ = "patch" + Foam::name(patchi);
+        bentry.type_ = "patch";
+        bentry.size_ = 0;
+        bentry.faces_ = labelList(one{}, patchi);
+    }
+}
+
+
 void Foam::PDRblock::adjustSizes()
 {
     // Adjust i-j-k addressing
@@ -438,9 +454,7 @@ void Foam::PDRblock::readBoundary(const dictionary& dict)
 
         // Save information for later access during mesh creation.
 
-        patches_.set(patchi, new boundaryEntry());
-
-        boundaryEntry& bentry = patches_[patchi];
+        boundaryEntry& bentry = patches_.emplace_set(patchi);
 
         bentry.name_ = patchName;
         bentry.type_ = patchType;
@@ -462,9 +476,7 @@ void Foam::PDRblock::readBoundary(const dictionary& dict)
 
     if (missed.size())
     {
-        patches_.append(new boundaryEntry());
-
-        boundaryEntry& bentry = patches_.last();
+        boundaryEntry& bentry = patches_.emplace_back();
 
         bentry.name_ = "defaultFaces";
         bentry.type_ = emptyPolyPatch::typeName;
@@ -523,24 +535,20 @@ Foam::PDRblock::PDRblock
 :
     PDRblock(dictionary::null, false)
 {
-    // Default boundaries with patchi == shapeFacei
-    patches_.resize(6);
-    for (label patchi=0; patchi < 6; ++patchi)
-    {
-        patches_.set(patchi, new boundaryEntry());
-
-        boundaryEntry& bentry = patches_[patchi];
-
-        bentry.name_ = "patch" + Foam::name(patchi);
-        bentry.type_ = "patch";
-        bentry.size_ = 0;
-        bentry.faces_ = labelList(one{}, patchi);
-    }
-
+    addDefaultPatches();
     reset(xgrid, ygrid, zgrid);
 }
 
 
+Foam::PDRblock::PDRblock(const boundBox& box, const labelVector& nCells)
+:
+    PDRblock(dictionary::null, false)
+{
+    addDefaultPatches();
+    reset(box, nCells);
+}
+
+
 Foam::PDRblock::PDRblock(const dictionary& dict, bool verboseOutput)
 :
     ijkMesh(),
@@ -631,6 +639,28 @@ void Foam::PDRblock::reset
 }
 
 
+void Foam::PDRblock::reset(const boundBox& box, const labelVector& nCells)
+{
+    grid_.x().reset(box.min().x(), box.max().x(), nCells.x());
+    grid_.y().reset(box.min().y(), box.max().y(), nCells.y());
+    grid_.z().reset(box.min().z(), box.max().z(), nCells.z());
+
+    adjustSizes();
+
+    // Adjust boundaries
+    for (boundaryEntry& bentry : patches_)
+    {
+        bentry.size_ = 0;
+
+        // Count patch faces
+        for (const label shapeFacei : bentry.faces_)
+        {
+            bentry.size_ += nBoundaryFaces(shapeFacei);
+        }
+    }
+}
+
+
 bool Foam::PDRblock::findCell(const point& pt, labelVector& pos) const
 {
     // Out-of-bounds is handled explicitly, for efficiency and consistency,

src/mesh/blockMesh/PDRblockMesh/PDRblock.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,8 +65,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef PDRblock_H
-#define PDRblock_H
+#ifndef Foam_PDRblock_H
+#define Foam_PDRblock_H
 
 #include "ijkMesh.H"
 #include "boundBox.H"
@@ -122,6 +122,9 @@ public:
         {
         public:
 
+            //- Reset with min/max range and number of divisions
+            void reset(const scalar low, const scalar upp, const label nCells);
+
             //- The location list is valid if it contains 2 or more points
             inline bool valid() const;
 
@@ -338,6 +341,9 @@ private:
             const UList<scalar>& pts
         );
 
+        //- Add zero-sized patches with default naming/types
+        void addDefaultPatches();
+
         //- Adjust sizing for updated grid points
         void adjustSizes();
 
@@ -440,6 +446,9 @@ public:
             const UList<scalar>& zgrid
         );
 
+        //- Construct from cube with specified griding
+        PDRblock(const boundBox& box, const labelVector& nCells);
+
         //- Construct from dictionary
         explicit PDRblock(const dictionary& dict, bool verboseOutput=false);
 
@@ -457,11 +466,14 @@ public:
             const UList<scalar>& zgrid
         );
 
+        //- Reset cube and mesh i-j-k sizing
+        void reset(const boundBox& box, const labelVector& nCells);
+
 
     // Access
 
         //- The grid point locations in the i,j,k (x,y,z) directions.
-        inline const Vector<location>& grid() const;
+        const Vector<location>& grid() const noexcept { return grid_; }
 
         //- Equivalent edge grading descriptors in (x,y,z) directions.
         Vector<gradingDescriptors> grading() const;
@@ -476,10 +488,10 @@ public:
         using ijkMesh::sizes;
 
         //- The mesh bounding box
-        inline const boundBox& bounds() const;
+        const boundBox& bounds() const noexcept { return bounds_; }
 
         //- The min/max edge length
-        inline const scalarMinMax& edgeLimits() const;
+        const scalarMinMax& edgeLimits() const noexcept { return edgeLimits_; }
 
         //- Cell size in x-direction at i position.
         inline scalar dx(const label i) const;

src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C

@@ -252,7 +252,7 @@ Foam::Ostream& Foam::PDRblock::blockMeshDict
                 *dummyTimePtr,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false  // no register
+                IOobject::NO_REGISTER
             )
         );
 
@@ -375,7 +375,7 @@ Foam::Ostream& Foam::PDRblock::blockMeshDict
                 *dummyTimePtr,
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false  // do not register
+                IOobject::NO_REGISTER
             );
 
             searchableSphere sphere(io, radialCentre, radialSizes);
@@ -796,7 +796,7 @@ void Foam::PDRblock::writeBlockMeshDict(const IOobject& io) const
             io.db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
-            false  // no register
+            IOobject::NO_REGISTER
         )
     );
 
@@ -831,7 +831,7 @@ Foam::PDRblock::createBlockMesh(const IOobject& io) const
             io.db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
-            false  // no register
+            IOobject::NO_REGISTER
         ),
         blockMeshDict()
     );

src/mesh/blockMesh/PDRblockMesh/PDRblockI.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,31 +48,31 @@ inline Foam::label Foam::PDRblock::location::nPoints() const
 
 inline bool Foam::PDRblock::location::contains(const scalar p) const
 {
-    return (scalarList::size() > 1 && first() <= p && p <= last());
+    return (scalarList::size() > 1 && front() <= p && p <= back());
 }
 
 
 inline const Foam::scalar& Foam::PDRblock::location::min() const
 {
-    return scalarList::empty() ? pTraits<scalar>::rootMax : first();
+    return scalarList::empty() ? pTraits<scalar>::rootMax : front();
 }
 
 
 inline const Foam::scalar& Foam::PDRblock::location::max() const
 {
-    return scalarList::empty() ? pTraits<scalar>::rootMin : last();
+    return scalarList::empty() ? pTraits<scalar>::rootMin : back();
 }
 
 
 inline Foam::scalar Foam::PDRblock::location::centre() const
 {
-    return scalarList::empty() ? 0 : (0.5*first() + 0.5*last());
+    return scalarList::empty() ? 0 : (0.5*front() + 0.5*back());
 }
 
 
 inline Foam::scalar Foam::PDRblock::location::length() const
 {
-    return scalarList::empty() ? 0 : mag(last() - first());
+    return scalarList::empty() ? 0 : mag(back() - front());
 }
 
 
@@ -107,12 +107,12 @@ inline Foam::scalar Foam::PDRblock::location::C(const label i) const
         #endif
 
         // "Halo" centre [-1] == x0 - 1/2 (x1 - x0)
-        return first() - 0.5*(operator[](1) - first());
+        return front() - 0.5*(operator[](1) - front());
     }
     else if (i > 1 && i == scalarList::size()-1)
     {
         // "Halo" centre [nCells] == xN + 1/2 (xN - xN1)
-        return last() + 0.5*(last() - operator[](scalarList::size()-2));
+        return back() + 0.5*(back() - operator[](scalarList::size()-2));
     }
 
     #ifdef FULLDEBUG
@@ -128,13 +128,13 @@ Foam::PDRblock::location::clip(const scalar& val) const
 {
     if (scalarList::size())
     {
-        if (val < first())
+        if (val < front())
         {
-            return first();
+            return front();
         }
-        else if (last() < val)
+        else if (back() < val)
         {
-            return last();
+            return back();
         }
     }
 
@@ -144,25 +144,6 @@ Foam::PDRblock::location::clip(const scalar& val) const
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline const Foam::Vector<Foam::PDRblock::location>&
-Foam::PDRblock::grid() const
-{
-    return grid_;
-}
-
-
-inline const Foam::scalarMinMax& Foam::PDRblock::edgeLimits() const
-{
-    return edgeLimits_;
-}
-
-
-inline const Foam::boundBox& Foam::PDRblock::bounds() const
-{
-    return bounds_;
-}
-
-
 inline Foam::scalar Foam::PDRblock::dx(const label i) const
 {
     return grid_.x().width(i);

src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -260,6 +260,29 @@ void Foam::PDRblock::gridControl::writeDict
 }
 
 
+void Foam::PDRblock::location::reset
+(
+    const scalar low,
+    const scalar upp,
+    const label nCells
+)
+{
+    scalarList& grid = *this;
+
+    grid.resize_nocopy(nCells+1);
+
+    grid.front() = low;
+    grid.back() = upp;
+
+    const scalar span = (upp - low);
+
+    for (label pointi = 1; pointi < nCells; ++pointi)
+    {
+        grid[pointi] = low + (pointi*span)/nCells;
+    }
+}
+
+
 Foam::scalarMinMax Foam::PDRblock::location::edgeLimits() const
 {
     scalarMinMax limits;

src/meshTools/coordSet/writers/vtk/foamVtkCoordSetWriter.H

@@ -201,8 +201,9 @@ public:
         void piece(const UPtrList<const pointField>& points);
 
 
-        //- Write processor ids for each line as CellData
-        //- (no-op in serial)
+        //- Write processor ids for each line as CellData or for each point
+        //- as PointData, depending on isPointData() state. No-op in serial.
+        //  Not implemented.
         bool writeProcIDs();
 
 

src/meshTools/output/vtk/patch/foamVtkGenericPatchWriter.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -159,6 +159,10 @@ public:
         //- Write processor ids for each poly as CellData
         bool writeProcIDs()
         {
+            if (this->isPointData())
+            {
+                return vtk::fileWriter::writeProcIDs(nLocalPoints_);
+            }
             return vtk::polyWriter::writeProcIDs(nLocalPolys_);
         }
 };

src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -585,6 +585,10 @@ void Foam::vtk::patchMeshWriter::writePatchIDs()
 
 bool Foam::vtk::patchMeshWriter::writeProcIDs()
 {
+    if (this->isPointData())
+    {
+        return vtk::fileWriter::writeProcIDs(nLocalPoints_);
+    }
     return vtk::fileWriter::writeProcIDs(nLocalPolys_);
 }
 

src/meshTools/output/vtk/patch/foamVtkPatchMeshWriter.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -209,8 +209,8 @@ public:
         //  Must be called within the CELL_DATA state.
         void writePatchIDs();
 
-        //- Write processor ids as CellData. This is no-op in serial.
-        //  Must be called within the CELL_DATA state.
+        //- Write processor ids for each line as CellData or for each point
+        //- as PointData, depending on isPointData() state. No-op in serial.
         bool writeProcIDs();
 
         //- Write processor neighbour ids as CellData. This is no-op in serial.

tutorials/lagrangian/kinematicParcelFoam/drippingChair/Allrun

@@ -3,19 +3,30 @@ cd "${0%/*}" || exit                                # Run from this directory
 . ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
 #------------------------------------------------------------------------------
 
+unset decompDict
+# decompDict=system/decomposeParDict-6
+# decompDict=system/decomposeParDict-7
+
 runApplication blockMesh
 
 runApplication makeFaMesh
 
-runApplication decomposePar
+runApplication -decompose-dict="$decompDict" decomposePar
 
-# For the processor partitioning
-runParallel -s finiteVolume foamToVTK -name VTK-parallel -time 0 \
-    -no-finite-area -no-internal -no-lagrangian -no-fields -with-ids
+# For ids and processor partitioning
+runParallel -s volume -decompose-dict="$decompDict" \
+    foamToVTK -name VTK-parallel -time 0 \
+    -no-finite-area -no-internal -no-lagrangian -no-fields -with-ids \
+    -patches filmWalls
 
-runParallel $(getApplication)
+runParallel -decompose-dict="$decompDict" $(getApplication)
 
-runParallel -s finiteArea foamToVTK -name VTK-parallel \
-    -no-boundary -no-internal -no-lagrangian
+if false
+then
+    # Not usually needed - uses areaWrite
+    runParallel -s area -decompose-dict="$decompDict" \
+        foamToVTK -name VTK-parallel \
+        -no-boundary -no-internal -no-lagrangian -with-ids
+fi
 
 #------------------------------------------------------------------------------

tutorials/lagrangian/kinematicParcelFoam/drippingChair/Allrun-serial

@@ -7,13 +7,18 @@ runApplication blockMesh
 
 runApplication makeFaMesh
 
-# For the cell ids etc
-runApplication -s finiteVolume.serial foamToVTK -name VTK-serial -time 0 \
-    -no-finite-area -no-internal -no-lagrangian -no-fields -with-ids
+# For ids and processor partitioning
+runApplication -s volume.serial foamToVTK -name VTK-serial -time 0 \
+    -no-finite-area -no-internal -no-lagrangian -no-fields -with-ids \
+    -patches filmWalls
 
 runApplication $(getApplication)
 
-runApplication -s finiteArea.serial foamToVTK -name VTK-serial \
-    -no-boundary -no-internal -no-lagrangian
+if false
+then
+    # Not usually needed - uses areaWrite
+    runApplication -s area.serial foamToVTK -name VTK-serial \
+        -no-boundary -no-internal -no-lagrangian
+fi
 
 #------------------------------------------------------------------------------

tutorials/lagrangian/kinematicParcelFoam/drippingChair/system/areaWrite

@@ -0,0 +1,25 @@
+// -*- C++ -*-
+// Use the areaWrite function object
+
+areaWrite
+{
+    type    areaWrite;
+    libs    (utilityFunctionObjects);
+    log     true;
+
+    writeControl    writeTime;
+    writeInterval   1;
+
+    // Fields to output (words or regex)
+    fields  (Uf_film hf_film pf_film);
+
+    surfaceFormat ensight;
+
+    formatOptions
+    {
+        default { format binary; }
+    }
+}
+
+
+// ************************************************************************* //

tutorials/lagrangian/kinematicParcelFoam/drippingChair/system/controlDict

@@ -52,5 +52,10 @@ regionFaMaxCo   1;
 
 maxDeltaT       0.1;
 
+functions
+{
+    #include "areaWrite"
+}
+
 
 // ************************************************************************* //

tutorials/lagrangian/kinematicParcelFoam/drippingChair/system/decomposeParDict-6

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2306                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 6;
+
+method          scotch;
+
+// ************************************************************************* //

tutorials/lagrangian/kinematicParcelFoam/drippingChair/system/decomposeParDict-7

@@ -0,0 +1,21 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2306                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      decomposeParDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+numberOfSubdomains 7;
+
+method          scotch;
+
+// ************************************************************************* //