Improved point-cell and cell-point topology methods

Goal

The goal of this change is to reduce the runtime of the pointCells() and cellPoints() methods in the primitiveMesh class, which calculate the point-to-cell and cell-to-point addressing, respectively.

Proposed Change

The pointCells() currently calculates the topology using the private method calcPointCells(), which utilizes the method labels() in the cell class. This method introduces significant overhead and reduces performance. The cellPoints() method currently calls the pointCell() method and inverts the result.

The proposed solution implements the cellPoints() method directly using a more efficient algorithm that makes use of the bitset data class to mark points that have already been found. Then, the pointCells() method is implemented by calling cellPoints() and inverting the result. As can be seen below, our implementation improves performance in all cases, even if the user only wants the point-cell topology.

Results

The speedup values were calculated by running both the current version and the new version on structured cube meshes of varying sizes. It can be seen that the speedup remains constant among all the different mesh sizes.
In the testing, the meshes began without either one of the topologies present. Then 3 different cases were tested:

1. The user wants only point-cell topology
   • Speedup = 1.4x
2. The user wants only cell-point topology
   • Speedup = 2.4x
3. The user wants both topologies
   • Speedup = 1.7x

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018 OpenCFD Ltd.
+    Copyright (C) 2018-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ License
 #include "primitiveMesh.H"
 #include "cell.H"
 #include "bitSet.H"
+#include "DynamicList.H"
 #include "ListOps.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -57,12 +58,51 @@ void Foam::primitiveMesh::calcCellPoints() const
             << "cellPoints already calculated"
             << abort(FatalError);
     }
-    else
+    else if (hasPointCells())
     {
         // Invert pointCells
         cpPtr_ = new labelListList(nCells());
         invertManyToMany(nCells(), pointCells(), *cpPtr_);
     }
+    else
+    {
+        // Calculate cell-point topology
+
+        cpPtr_ = new labelListList(nCells());
+        auto& cellPointAddr = *cpPtr_;
+
+        const cellList& cellLst = cells();
+        const faceList& faceLst = faces();
+
+        // Tracking (only use each point id once)
+        bitSet usedPoints(nPoints());
+
+        // Vertex labels for the current cell
+        DynamicList<label> currPoints(256);
+
+        const label loopLen = nCells();
+
+        for (label celli = 0; celli < loopLen; ++celli)
+        {
+            // Clear any previous contents
+            usedPoints.unset(currPoints);
+            currPoints.clear();
+
+            for (const label facei : cellLst[celli])
+            {
+                for (const label pointi : faceLst[facei])
+                {
+                    // Only once for each point id
+                    if (usedPoints.set(pointi))
+                    {
+                        currPoints.push_back(pointi);
+                    }
+                }
+            }
+
+            cellPointAddr[celli] = currPoints;  // NB: unsorted
+        }
+    }
 }