Merge branch 'primitiveMeshOptimization' into 'develop'

ENH: improved point-cell and cell-point topology methods (#2715)

See merge request Development/openfoam!597

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
+        }
+    }
 }
 
 

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,128 @@ void Foam::primitiveMesh::calcPointCells() const
             << "pointCells already calculated"
             << abort(FatalError);
     }
-    else
+    else if (hasCellPoints())
+    {
+        // Invert cellPoints
+        pcPtr_ = new labelListList(nPoints());
+        invertManyToMany(nPoints(), cellPoints(), *pcPtr_);
+    }
+    else if (hasPointFaces())
     {
-        const cellList& cf = cells();
+        // Calculate point-cell from point-face information
+
+        const labelList& own = faceOwner();
+        const labelList& nei = faceNeighbour();
+        const labelListList& pFaces = pointFaces();
 
-        // Count number of cells per point
+        // Tracking (only use each cell id once)
+        bitSet usedCells(nCells());
 
-        labelList npc(nPoints(), Zero);
+        // Cell ids for the point currently being processed
+        DynamicList<label> currCells(256);
 
-        forAll(cf, celli)
+        const label loopLen = nPoints();
+
+        pcPtr_ = new labelListList(nPoints());
+        auto& pointCellAddr = *pcPtr_;
+
+        for (label pointi = 0; pointi < loopLen; ++pointi)
         {
-            const labelList curPoints = cf[celli].labels(faces());
+            // Clear any previous contents
+            usedCells.unset(currCells);
+            currCells.clear();
 
-            forAll(curPoints, pointi)
+            for (const label facei : pFaces[pointi])
             {
-                label ptI = curPoints[pointi];
+                // Owner cell - only allow one occurance
+                if (usedCells.set(own[facei]))
+                {
+                    currCells.push_back(own[facei]);
+                }
+
+                // Neighbour cell - only allow one occurance
+                if (facei < nInternalFaces())
+                {
+                    if (usedCells.set(nei[facei]))
+                    {
+                        currCells.push_back(nei[facei]);
+                    }
+                }
+            }
+
+            pointCellAddr[pointi] = currCells;  // NB: unsorted
+        }
+    }
+    else
+    {
+        // Calculate point-cell topology
+
+        const cellList& cellLst = cells();
+        const faceList& faceLst = faces();
+
+        // Tracking (only use each point id once)
+        bitSet usedPoints(nPoints());
+
+        // Which of usedPoints needs to be unset [faster]
+        DynamicList<label> currPoints(256);
+
+        const label loopLen = nCells();
+
+        // Step 1: count number of cells per point
+
+        labelList pointCount(nPoints(), Zero);
 
-                npc[ptI]++;
+        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);  // Needed for cleanup
+                        ++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
+            usedPoints.unset(currPoints);
+            currPoints.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 (usedPoints.set(pointi))
+                    {
+                        currPoints.push_back(pointi);  // Needed for cleanup
+                        pointCellAddr[pointi][pointCount[pointi]++] = celli;
+                    }
+                }
             }
         }
     }
@@ -154,12 +234,8 @@ const Foam::labelList& Foam::primitiveMesh::pointCells
         if (storage.size() > 1)
         {
             std::sort(storage.begin(), storage.end());
-
             auto last = std::unique(storage.begin(), storage.end());
-
-            const label newLen = label(last - storage.begin());
-
-            storage.resize(newLen);
+            storage.resize(label(last - storage.begin()));
         }
 
         return storage;