ENH: improve efficiency of point-cell/cell-point construction (#2715)

- calculate point-cell from point-face if available

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,49 @@ 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();
+
+        // Vertex labels for the current cell
+        DynamicList<label> vertices(256);
+
+
+        const label loopLen = 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;
+        }
+    }
 }
 
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,14 +29,13 @@ License
 #include "primitiveMesh.H"
 #include "cell.H"
 #include "bitSet.H"
+#include "DynamicList.H"
 #include "ListOps.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::primitiveMesh::calcPointCells() const
 {
-    // Loop through cells and mark up points
-
     if (debug)
     {
         Pout<< "primitiveMesh::calcPointCells() : "
@@ -59,48 +59,121 @@ void Foam::primitiveMesh::calcPointCells() const
             << "pointCells already calculated"
             << abort(FatalError);
     }
+    else if (hasCellPoints())
+    {
+        // Invert cellPoints
+        pcPtr_ = new labelListList(nPoints());
+        invertManyToMany(nPoints(), cellPoints(), *pcPtr_);
+    }
+    else if (hasPointFaces())
+    {
+        // Calculate point-cell from point-face information
+
+        const labelList& own = faceOwner();
+        const labelList& nei = faceNeighbour();
+        const labelListList& pFaces = pointFaces();
+
+        const label loopLen = nPoints();
+
+        pcPtr_ = new labelListList(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 < 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;
+        }
+    }
     else
     {
-        const cellList& cf = cells();
+        // Calculate point-cell topology
 
-        // Count number of cells per point
+        const cellList& cellLst = cells();
+        const faceList& faceLst = faces();
 
-        labelList npc(nPoints(), Zero);
+        // Vertex labels for the current cell
+        // Note: since the number of points per cell is relatively small,
+        // using a bitSet to track does not help much.
+        // Using a labelHashSet to track has too many allocations.
+        DynamicList<label> vertices(256);
 
-        forAll(cf, celli)
+
+        const label loopLen = nCells();
+
+        // Step 1: count number of cells per point
+
+        labelList pointCount(nPoints(), Zero);
+
+        for (label celli = 0; celli < loopLen; ++celli)
         {
-            const labelList curPoints = cf[celli].labels(faces());
+            // Clear any previous contents
+            vertices.clear();
 
-            forAll(curPoints, pointi)
+            for (const label facei : cellLst[celli])
             {
-                label ptI = curPoints[pointi];
-
-                npc[ptI]++;
+                for (const label pointi : faceLst[facei])
+                {
+                    // Only once for each point id
+                    if (!vertices.contains(pointi))
+                    {
+                        vertices.push_back(pointi);
+                        ++pointCount[pointi];
+                    }
+                }
             }
         }
 
 
-        // Size and fill cells per point
+        // Step 2: set sizing, reset counters
 
-        pcPtr_ = new labelListList(npc.size());
-        labelListList& pointCellAddr = *pcPtr_;
+        pcPtr_ = new labelListList(nPoints());
+        auto& pointCellAddr = *pcPtr_;
 
         forAll(pointCellAddr, pointi)
         {
-            pointCellAddr[pointi].setSize(npc[pointi]);
+            pointCellAddr[pointi].resize_nocopy(pointCount[pointi]);
+            pointCount[pointi] = 0;
         }
-        npc = 0;
 
 
-        forAll(cf, celli)
+        // Step 3: fill in values. Logic as per step 1
+        for (label celli = 0; celli < loopLen; ++celli)
         {
-            const labelList curPoints = cf[celli].labels(faces());
+            // Clear any previous contents
+            vertices.clear();
 
-            forAll(curPoints, pointi)
+            for (const label facei : cellLst[celli])
             {
-                label ptI = curPoints[pointi];
-
-                pointCellAddr[ptI][npc[ptI]++] = 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;
+                    }
+                }
             }
         }
     }