From 709472a7bbe49ddfe7e522f0153a599f00cc8e4e Mon Sep 17 00:00:00 2001
From: alonzameret <alonzameret@gmail.com>
Date: Mon, 27 Feb 2023 15:32:32 +0200
Subject: [PATCH] ENH: improve efficiency of point-cell/cell-point construction
 (#2715)

- reworked simpler looping, reinstate pointCells() calculation
  but with markup as per cellPoints()
---
 .../primitiveMesh/primitiveMeshCellPoints.C   | 61 ++++++---------
 .../primitiveMesh/primitiveMeshPointCells.C   | 77 ++++++++++++++++++-
 2 files changed, 100 insertions(+), 38 deletions(-)

diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C
index 9748e8f9a1..a0e1f3d7e6 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C
@@ -29,6 +29,7 @@ License
 #include "primitiveMesh.H"
 #include "cell.H"
 #include "bitSet.H"
+#include "DynamicList.H"
 #include "ListOps.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -67,52 +68,40 @@ void Foam::primitiveMesh::calcCellPoints() const
     {
         // Calculate cell-point topology
 
-        const cellList& cL = cells();
-        const faceList& fL = faces();
-
         cpPtr_ = new labelListList(nCells());
-        labelListList& cellPoints = *cpPtr_;
+        auto& cellPointAddr = *cpPtr_;
 
-        // Mark points we have found as to not count them twice for the same cell
-        Foam::bitSet foundPoint(nPoints(), false);
+        const cellList& cellLst = cells();
+        const faceList& faceLst = faces();
 
-        forAll (cL, cellIdx)
+        // For tracking (only use each point id once)
+        bitSet usedPoints(nPoints());
+
+        // Vertex labels for the current cell
+        DynamicList<label> vertices(256);
+
+
+        const label loopLen = nCells();
+
+        for (label celli = 0; celli < loopLen; ++celli)
         {
-            // List of faces on the cell
-            const Foam::labelList& cellFaces = cL[cellIdx];
-            // List that will contain the cell's points
-            Foam::labelList& neiPoints = cellPoints[cellIdx];
+            // Clear any previous contents
+            usedPoints.unset(vertices);
+            vertices.clear();
 
-            // Over-Estimate number of neighbour points
-            int pointCnt = 0;
-            for (label faceIdx : cellFaces)
+            for (const label facei : cellLst[celli])
             {
-                pointCnt += fL[faceIdx].size();
-            }
-            neiPoints.resize(pointCnt);
-
-            // Find the neighbour points
-            pointCnt = 0;
-            for (label faceIdx : cellFaces)
-            {
-                for (label pointIdx : fL[faceIdx])
+                for (const label pointi : faceLst[facei])
                 {
-                    // This check guarantees every point is added exactly once
-                    if (foundPoint[pointIdx]) continue;
-
-                    foundPoint[pointIdx] = true;
-                    neiPoints[pointCnt++] = pointIdx;
+                    // Only once for each point id
+                    if (usedPoints.set(pointi))
+                    {
+                        vertices.push_back(pointi);
+                    }
                 }
             }
 
-            // Resize to actual value
-            neiPoints.resize(pointCnt);
-
-            // Clear bitSet before next iteration
-            for (label pointIdx : neiPoints)
-            {
-                foundPoint[pointIdx] = false;
-            }
+            cellPointAddr[celli] = vertices;  // unsorted
         }
     }
 }
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C
index ed60da9e59..9bd875e18d 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C
@@ -29,6 +29,7 @@ License
 #include "primitiveMesh.H"
 #include "cell.H"
 #include "bitSet.H"
+#include "DynamicList.H"
 #include "ListOps.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -58,13 +59,85 @@ void Foam::primitiveMesh::calcPointCells() const
             << "pointCells already calculated"
             << abort(FatalError);
     }
-    else
+    else if (hasCellPoints())
     {
         // Invert cellPoints
-        (void)cellPoints();
         pcPtr_ = new labelListList(nPoints());
         invertManyToMany(nPoints(), cellPoints(), *pcPtr_);
     }
+    else
+    {
+        // Calculate point-cell topology
+
+        const cellList& cellLst = cells();
+        const faceList& faceLst = faces();
+
+        // For tracking (only use each point id once)
+        bitSet usedPoints(nPoints());
+
+        // Vertex labels for the current cell
+        DynamicList<label> vertices(256);
+
+        const label loopLen = nCells();
+
+        // Step 1: count number of cells per point
+
+        labelList pointCount(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
+
+        pcPtr_ = new labelListList(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;
+                    }
+                }
+            }
+        }
+    }
 }