songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

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

Browse Source 
- for OpenFOAM-v2212 and earlier cellPoints() were constructed
  from pointCells(), but this is slower than constructing
  pointCells() from cellPoints().

  Some of the slowness is due to allocations associated with
  cells::labels(), but a large amount of slowness is the duplicate
  point avoidance. Since this is being done for many points/cells
  at once, using a bitSet for managing the duplicates amortizes
  quickly

- now construct cellPoints() from cached pointCells(), otherwise
  construct manually (using bitSet/DynamicList for bookkeeping)

- construct pointCells() from cached cellPoints(), or cached
  pointFaces(), otherwise manually.

Code Contribution: Alon Zameret

Co-authored-by: Mark Olesen
This commit is contained in:
alonzameret 2023-02-27 15:32:32 +02:00 committed by Mark Olesen
parent 820f93c3b4
commit 61d610574c
2 changed files with 145 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellPoints.C
View File

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

130
src/OpenFOAM/meshes/primitiveMesh/primitiveMeshPointCells.C
View File

@ -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 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();
// Tracking (only use each cell id once)
bitSet usedCells(nCells());
// Cell ids for the point currently being processed
DynamicList<label> currCells(256);
const label loopLen = nPoints();
pcPtr_ = new labelListList(nPoints());
auto& pointCellAddr = *pcPtr_;
for (label pointi = 0; pointi < loopLen; ++pointi)
{
// Clear any previous contents
usedCells.unset(currCells);
currCells.clear();
for (const label facei : pFaces[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
{
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);
// Tracking (only use each point id once)
bitSet usedPoints(nPoints());
forAll(cf, celli)
// 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);
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];
npc[ptI]++;
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;