From 820f93c3b4cde81fe4c99366577bc4e655b38878 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 3 Mar 2023 12:27:44 +0100
Subject: [PATCH] ENH: simplify box mesh construction (PDRblockMesh)

ENH: construct or reset boundBox from zero_one for a 0-1 unit cube

ENH: add some testing timings (for !596 and #2715)
---
 applications/test/boundBox2/Test-boundBox2.C  |   2 +-
 .../test/polyMeshGeom-speed1/Make/files       |   3 +
 .../test/polyMeshGeom-speed1/Make/options     |   5 +
 .../Test-polyMeshGeom-speed1.C                | 804 ++++++++++++++++++
 .../Test-triangleIntersection.C               |   3 +-
 src/OpenFOAM/meshes/boundBox/boundBox.H       |   8 +-
 src/OpenFOAM/meshes/boundBox/boundBoxI.H      |  16 +-
 src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H   |   4 +-
 src/OpenFOAM/meshes/ijkMesh/ijkMesh.H         |   4 +-
 .../meshes/treeBoundBox/treeBoundBox.H        |   5 +-
 .../meshes/treeBoundBox/treeBoundBoxI.H       |   8 +-
 src/mesh/blockMesh/PDRblockMesh/PDRblock.C    |  72 +-
 src/mesh/blockMesh/PDRblockMesh/PDRblock.H    |  24 +-
 .../PDRblockMesh/PDRblockBlockMesh.C          |   8 +-
 src/mesh/blockMesh/PDRblockMesh/PDRblockI.H   |  43 +-
 .../blockMesh/PDRblockMesh/PDRblockLocation.C |  25 +-
 16 files changed, 960 insertions(+), 74 deletions(-)
 create mode 100644 applications/test/polyMeshGeom-speed1/Make/files
 create mode 100644 applications/test/polyMeshGeom-speed1/Make/options
 create mode 100644 applications/test/polyMeshGeom-speed1/Test-polyMeshGeom-speed1.C

diff --git a/applications/test/boundBox2/Test-boundBox2.C b/applications/test/boundBox2/Test-boundBox2.C
index 3c88e70d51..4d929dd7a4 100644
--- a/applications/test/boundBox2/Test-boundBox2.C
+++ b/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
diff --git a/applications/test/polyMeshGeom-speed1/Make/files b/applications/test/polyMeshGeom-speed1/Make/files
new file mode 100644
index 0000000000..e6c64cc479
--- /dev/null
+++ b/applications/test/polyMeshGeom-speed1/Make/files
@@ -0,0 +1,3 @@
+Test-polyMeshGeom-speed1.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-polyMeshGeom-speed1
diff --git a/applications/test/polyMeshGeom-speed1/Make/options b/applications/test/polyMeshGeom-speed1/Make/options
new file mode 100644
index 0000000000..b4f99e5295
--- /dev/null
+++ b/applications/test/polyMeshGeom-speed1/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+    -I$(LIB_SRC)/mesh/blockMesh/lnInclude
+
+EXE_LIBS = \
+    -lblockMesh
diff --git a/applications/test/polyMeshGeom-speed1/Test-polyMeshGeom-speed1.C b/applications/test/polyMeshGeom-speed1/Test-polyMeshGeom-speed1.C
new file mode 100644
index 0000000000..4c268797dd
--- /dev/null
+++ b/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;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/test/triangleIntersection/Test-triangleIntersection.C b/applications/test/triangleIntersection/Test-triangleIntersection.C
index 1aca11601c..6573ec1b3e 100644
--- a/applications/test/triangleIntersection/Test-triangleIntersection.C
+++ b/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);
diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H
index 1373c5d754..554c93285c 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBox.H
+++ b/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);
 
diff --git a/src/OpenFOAM/meshes/boundBox/boundBoxI.H b/src/OpenFOAM/meshes/boundBox/boundBoxI.H
index 94bd1ca6f7..4b900ae1c5 100644
--- a/src/OpenFOAM/meshes/boundBox/boundBoxI.H
+++ b/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;
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H b/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H
index dacef4869a..1409d3b967 100644
--- a/src/OpenFOAM/meshes/ijkMesh/ijkAddressing.H
+++ b/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"
diff --git a/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H b/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H
index 8da2879a20..dc73b1fee0 100644
--- a/src/OpenFOAM/meshes/ijkMesh/ijkMesh.H
+++ b/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"
 
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
index 5b7315b3bf..a30a872250 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H
+++ b/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);
 
diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H
index 1db71c74bd..08bff3acfa 100644
--- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBoxI.H
+++ b/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)
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C b/src/mesh/blockMesh/PDRblockMesh/PDRblock.C
index c1b6f3b6f0..5e27e50075 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblock.C
+++ b/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,
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H b/src/mesh/blockMesh/PDRblockMesh/PDRblock.H
index d9cffa0131..a9096929f8 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblock.H
+++ b/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;
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C
index 5d8032a66e..13f4d8bba5 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblockBlockMesh.C
+++ b/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()
     );
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H b/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H
index d05c279dd6..22dc919fe1 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblockI.H
+++ b/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);
diff --git a/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C b/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C
index 05138a4c2a..b5eb724462 100644
--- a/src/mesh/blockMesh/PDRblockMesh/PDRblockLocation.C
+++ b/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;