Issue 3215 point connected cells

This commit is contained in:
Mattijs Janssens
2024-09-18 11:07:41 +00:00
committed by Andrew Heather
parent 1d6396dd3f
commit d72f51ac23
7 changed files with 495 additions and 201 deletions

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020,2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,60 +27,137 @@ License
\*---------------------------------------------------------------------------*/
#include "nearWallDist.H"
#include "fvMesh.H"
#include "cellDistFuncs.H"
#include "wallFvPatch.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::nearWallDist::calculate()
{
cellDistFuncs wallUtils(mesh_);
// Get patch ids of walls
labelHashSet wallPatchIDs(wallUtils.getPatchIDs<wallPolyPatch>());
// Size neighbours array for maximum possible
DynamicList<label> neighbours(wallUtils.maxPatchSize(wallPatchIDs));
// Correct all cells with face on wall
const cellDistFuncs wallUtils(mesh_);
const volVectorField& cellCentres = mesh_.C();
forAll(mesh_.boundary(), patchi)
if (cellDistFuncs::useCombinedWallPatch)
{
fvPatchScalarField& ypatch = operator[](patchi);
// Collect indices of wall patches
const fvPatch& patch = mesh_.boundary()[patchi];
DynamicList<label> wallPatchIDs(mesh_.boundary().size());
label nWalls = 0;
if (isA<wallFvPatch>(patch))
forAll(mesh_.boundary(), patchi)
{
const polyPatch& pPatch = patch.patch();
if (isA<wallFvPatch>(mesh_.boundary()[patchi]))
{
wallPatchIDs.append(patchi);
nWalls += mesh_.boundary()[patchi].size();
}
else
{
// Set distance to 0
operator[](patchi) = 0.0;
}
}
const labelUList& faceCells = patch.faceCells();
// Check cells with face on wall
// Collect all mesh faces of wall patches
DynamicList<label> faceLabels(nWalls);
for (const label patchi : wallPatchIDs)
{
const fvPatch& patch = mesh_.boundary()[patchi];
const auto& pp = patch.patch();
forAll(patch, i)
{
faceLabels.append(pp.start()+i);
}
}
const uindirectPrimitivePatch wallPatch
(
UIndirectList<face>(mesh_.faces(), faceLabels),
mesh_.points()
);
DynamicList<label> neighbours;
nWalls = 0;
for (const label patchi : wallPatchIDs)
{
const fvPatch& patch = mesh_.boundary()[patchi];
const labelUList& faceCells = patch.patch().faceCells();
fvPatchScalarField& ypatch = operator[](patchi);
forAll(patch, patchFacei)
{
wallUtils.getPointNeighbours(pPatch, patchFacei, neighbours);
// Get point connected neighbours (in wallPatch indices!)
wallUtils.getPointNeighbours(wallPatch, nWalls, neighbours);
label minFacei = -1;
ypatch[patchFacei] = wallUtils.smallestDist
(
cellCentres[faceCells[patchFacei]],
pPatch,
wallPatch,
neighbours,
minFacei
);
nWalls++;
}
}
else
}
else
{
// Get patch ids of walls
const labelHashSet wallPatchIDs(wallUtils.getPatchIDs<wallPolyPatch>());
// Size neighbours array for maximum possible
DynamicList<label> neighbours(wallUtils.maxPatchSize(wallPatchIDs));
// Correct all cells with face on wall
forAll(mesh_.boundary(), patchi)
{
ypatch = 0.0;
fvPatchScalarField& ypatch = operator[](patchi);
const fvPatch& patch = mesh_.boundary()[patchi];
if (isA<wallFvPatch>(patch))
{
const polyPatch& pPatch = patch.patch();
const labelUList& faceCells = patch.faceCells();
// Check cells with face on wall
forAll(patch, patchFacei)
{
wallUtils.getPointNeighbours
(
pPatch,
patchFacei,
neighbours
);
label minFacei = -1;
ypatch[patchFacei] = wallUtils.smallestDist
(
cellCentres[faceCells[patchFacei]],
pPatch,
neighbours,
minFacei
);
}
}
else
{
ypatch = 0.0;
}
}
}
}

View File

@ -234,19 +234,34 @@ void Foam::wallDistAddressing::correct(volScalarField& y)
if (correctWalls_)
{
cellToWallFace.reserve(nWalls);
correctBoundaryFaceCells
(
patchSet,
y,
cellToWallFace
);
correctBoundaryPointCells
(
patchSet,
y,
cellToWallFace
);
if (cellDistFuncs::useCombinedWallPatch)
{
// Correct across multiple patches
correctBoundaryCells
(
patchIDs_,
true, // do point-connected cells as well
y,
cellToWallFace
);
}
else
{
// Optional backwards compatibility
correctBoundaryFaceCells
(
patchSet,
y,
cellToWallFace
);
correctBoundaryPointCells
(
patchSet,
y,
cellToWallFace
);
}
}
// Make sure boundary values are up to date

View File

@ -29,6 +29,8 @@ License
#include "cellDistFuncs.H"
#include "polyMesh.H"
#include "polyBoundaryMesh.H"
#include "uindirectPrimitivePatch.H"
#include "registerSwitch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -37,6 +39,16 @@ namespace Foam
defineTypeNameAndDebug(cellDistFuncs, 0);
}
bool Foam::cellDistFuncs::useCombinedWallPatch = true;
registerInfoSwitch
(
"useCombinedWallPatch",
bool,
Foam::cellDistFuncs::useCombinedWallPatch
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::cellDistFuncs::cellDistFuncs(const polyMesh& mesh)
@ -56,134 +68,6 @@ Foam::labelHashSet Foam::cellDistFuncs::getPatchIDs
}
// Return smallest true distance from p to any of wallFaces.
// Note that even if normal hits face we still check other faces.
// Note that wallFaces is untruncated and we explicitly pass in size.
Foam::scalar Foam::cellDistFuncs::smallestDist
(
const point& p,
const polyPatch& patch,
const labelUList& wallFaces,
label& minFacei
) const
{
const pointField& points = patch.points();
scalar minDist = GREAT;
minFacei = -1;
for (const label patchFacei : wallFaces)
{
const pointHit curHit = patch[patchFacei].nearestPoint(p, points);
if (curHit.distance() < minDist)
{
minDist = curHit.distance();
minFacei = patch.start() + patchFacei;
}
}
return minDist;
}
// Get point neighbours of facei (including facei). Returns number of faces.
// Note: does not allocate storage but does use linear search to determine
// uniqueness. For polygonal faces this might be quite inefficient.
void Foam::cellDistFuncs::getPointNeighbours
(
const primitivePatch& patch,
const label patchFacei,
DynamicList<label>& neighbours
) const
{
neighbours.clear();
// Add myself
neighbours.append(patchFacei);
// Add all face neighbours
const labelList& faceNeighbours = patch.faceFaces()[patchFacei];
for (const label nbr : faceNeighbours)
{
neighbours.push_uniq(nbr);
}
// Add all point-only neighbours by linear searching in edge neighbours.
// Assumes that point-only neighbours are not using multiple points on
// face.
const face& f = patch.localFaces()[patchFacei];
forAll(f, fp)
{
label pointi = f[fp];
const labelList& pointNbs = patch.pointFaces()[pointi];
for (const label facei : pointNbs)
{
// Check for facei in edge-neighbours part of neighbours
neighbours.push_uniq(facei);
}
}
if (debug)
{
// Check for duplicates
// Use hashSet to determine nbs.
labelHashSet nbs(4*f.size());
forAll(f, fp)
{
const labelList& pointNbs = patch.pointFaces()[f[fp]];
nbs.insert(pointNbs);
}
// Subtract ours.
for (const label nb : neighbours)
{
if (!nbs.found(nb))
{
SeriousErrorInFunction
<< "getPointNeighbours : patchFacei:" << patchFacei
<< " verts:" << f << endl;
forAll(f, fp)
{
SeriousErrorInFunction
<< "point:" << f[fp] << " pointFaces:"
<< patch.pointFaces()[f[fp]] << endl;
}
for (const label facei : neighbours)
{
SeriousErrorInFunction
<< "fast nbr:" << facei
<< endl;
}
FatalErrorInFunction
<< "Problem: fast pointNeighbours routine included " << nb
<< " which is not in proper neighbour list " << nbs.toc()
<< abort(FatalError);
}
nbs.erase(nb);
}
if (nbs.size())
{
FatalErrorInFunction
<< "Problem: fast pointNeighbours routine did not find "
<< nbs.toc() << abort(FatalError);
}
}
}
// size of largest patch (out of supplied subset of patches)
Foam::label Foam::cellDistFuncs::maxPatchSize
(
@ -275,7 +159,7 @@ void Foam::cellDistFuncs::correctBoundaryFaceCells
);
// Store wallCell and its nearest neighbour
nearestFace.insert(celli, minFacei);
nearestFace.insert(celli, patch.start()+minFacei);
}
}
}
@ -348,7 +232,7 @@ void Foam::cellDistFuncs::correctBoundaryPointCells
);
// Store wallCell and its nearest neighbour
nearestFace.insert(celli, minFacei);
nearestFace.insert(celli, patch.start()+minFacei);
}
}
}
@ -357,4 +241,142 @@ void Foam::cellDistFuncs::correctBoundaryPointCells
}
void Foam::cellDistFuncs::correctBoundaryCells
(
const labelList& patchIDs,
const bool doPointCells,
scalarField& wallDistCorrected,
Map<label>& nearestFace
) const
{
label nWalls = 0;
{
for (const label patchi : patchIDs)
{
nWalls += mesh().boundaryMesh()[patchi].size();
}
}
DynamicList<label> faceLabels(nWalls);
{
for (const label patchi : patchIDs)
{
const auto& patch = mesh().boundaryMesh()[patchi];
forAll(patch, i)
{
faceLabels.append(patch.start()+i);
}
}
}
const uindirectPrimitivePatch wallPatch
(
UIndirectList<face>(mesh().faces(), faceLabels),
mesh_.points()
);
// Correct all cells with face on wall
const vectorField& cellCentres = mesh().cellCentres();
DynamicList<label> neighbours;
nWalls = 0;
for (const label patchi : patchIDs)
{
const auto& patch = mesh().boundaryMesh()[patchi];
const auto areaFraction(patch.areaFraction());
const labelUList& faceCells = patch.faceCells();
// Check cells with face on wall
forAll(patch, patchFacei)
{
if (areaFraction && (areaFraction()[patchFacei] <= 0.5))
{
// For cyclicACMI: more cyclic than wall
}
else
{
getPointNeighbours(wallPatch, nWalls, neighbours);
const label celli = faceCells[patchFacei];
label minFacei = -1;
wallDistCorrected[celli] = smallestDist
(
cellCentres[celli],
wallPatch,
neighbours,
minFacei
);
// Store wallCell and its nearest neighbour
nearestFace.insert(celli, nWalls+minFacei);
}
nWalls++;
}
}
// Correct all cells with a point on the wall
if (doPointCells)
{
const auto& meshPoints = wallPatch.meshPoints();
const auto& localFaces = wallPatch.localFaces();
bitSet isWallPoint(meshPoints.size(), true);
nWalls = 0;
for (const label patchi : patchIDs)
{
const auto& patch = mesh().boundaryMesh()[patchi];
const auto areaFraction(patch.areaFraction());
// Check cells with face on wall
forAll(patch, patchFacei)
{
if (areaFraction && (areaFraction()[patchFacei] <= 0.5))
{
// For cyclicACMI: more cyclic than wall
isWallPoint.unset(localFaces[nWalls]);
}
nWalls++;
}
}
const auto& pointFaces = wallPatch.pointFaces();
for (const label patchPointi : isWallPoint)
{
const label verti = meshPoints[patchPointi];
const labelList& neighbours = mesh().pointCells(verti);
for (const label celli : neighbours)
{
if (!nearestFace.found(celli))
{
const labelList& wallFaces = pointFaces[patchPointi];
label minFacei = -1;
wallDistCorrected[celli] = smallestDist
(
cellCentres[celli],
wallPatch,
wallFaces,
minFacei
);
// Store wallCell and its nearest neighbour
nearestFace.insert(celli, wallPatch.addressing()[minFacei]);
}
}
}
}
}
// ************************************************************************* //

View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020,2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -79,8 +79,16 @@ class cellDistFuncs
public:
// Static Data Members
ClassName("cellDistFuncs");
//- Use combined-wall-patches wall distance v.s. v2406 per-patch
//- distance. Default is true
static bool useCombinedWallPatch;
// Constructors
//- Construct from mesh
@ -103,22 +111,22 @@ public:
template<class Type>
labelHashSet getPatchIDs() const;
//- Calculate smallest true distance (and face index)
//- Calculate smallest true distance (and patch face index)
// from pt to faces wallFaces.
// For efficiency reasons we still pass in patch instead of extracting
// it from mesh_
template<class PatchType>
scalar smallestDist
(
const point& p,
const polyPatch& patch,
const PatchType& patch,
const labelUList& wallFaces,
label& meshFacei
label& patchFacei
) const;
//- Get faces sharing point with face on patch
template<class PatchType>
void getPointNeighbours
(
const primitivePatch&,
const PatchType&,
const label patchFacei,
DynamicList<label>&
) const;
@ -147,6 +155,17 @@ public:
scalarField& wallDistCorrected,
Map<label>& nearestFace
) const;
//- Correct all cells connected to any of the patches in patchIDs. Sets
// - cell values in wallDistCorrected
// - (mesh) face that contains the nearest point
void correctBoundaryCells
(
const labelList& patchIDs,
const bool doPointCells,
scalarField& wallDistCorrected,
Map<label>& nearestFace
) const;
};

View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2011-2016,2024 OpenFOAM Foundation
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -50,4 +50,135 @@ Foam::labelHashSet Foam::cellDistFuncs::getPatchIDs() const
}
template<class PatchType>
Foam::scalar Foam::cellDistFuncs::smallestDist
(
const point& p,
const PatchType& patch,
const labelUList& wallFaces,
label& minFacei
) const
{
// Return smallest true distance from p to any of wallFaces.
// Note that even if normal hits face we still check other faces.
const pointField& points = patch.points();
scalar minDist = GREAT;
minFacei = -1;
for (const label patchFacei : wallFaces)
{
const pointHit curHit = patch[patchFacei].nearestPoint(p, points);
if (curHit.distance() < minDist)
{
minDist = curHit.distance();
minFacei = patchFacei;
}
}
return minDist;
}
template<class PatchType>
void Foam::cellDistFuncs::getPointNeighbours
(
const PatchType& patch,
const label patchFacei,
DynamicList<label>& neighbours
) const
{
// Get point neighbours of facei (including facei). Returns number of faces.
// Note: does not allocate storage but does use linear search to determine
// uniqueness. For polygonal faces this might be quite inefficient.
neighbours.clear();
// Add myself
neighbours.append(patchFacei);
// Add all face neighbours
const labelList& faceNeighbours = patch.faceFaces()[patchFacei];
for (const label nbr : faceNeighbours)
{
neighbours.push_uniq(nbr);
}
// Add all point-only neighbours by linear searching in edge neighbours.
// Assumes that point-only neighbours are not using multiple points on
// face.
const face& f = patch.localFaces()[patchFacei];
forAll(f, fp)
{
label pointi = f[fp];
const labelList& pointNbs = patch.pointFaces()[pointi];
for (const label facei : pointNbs)
{
// Check for facei in edge-neighbours part of neighbours
neighbours.push_uniq(facei);
}
}
if (debug)
{
// Check for duplicates
// Use hashSet to determine nbs.
labelHashSet nbs(4*f.size());
forAll(f, fp)
{
const labelList& pointNbs = patch.pointFaces()[f[fp]];
nbs.insert(pointNbs);
}
// Subtract ours.
for (const label nb : neighbours)
{
if (!nbs.found(nb))
{
SeriousErrorInFunction
<< "getPointNeighbours : patchFacei:" << patchFacei
<< " verts:" << f << endl;
forAll(f, fp)
{
SeriousErrorInFunction
<< "point:" << f[fp] << " pointFaces:"
<< patch.pointFaces()[f[fp]] << endl;
}
for (const label facei : neighbours)
{
SeriousErrorInFunction
<< "fast nbr:" << facei
<< endl;
}
FatalErrorInFunction
<< "Problem: fast pointNeighbours routine included " << nb
<< " which is not in proper neighbour list " << nbs.toc()
<< abort(FatalError);
}
nbs.erase(nb);
}
if (nbs.size())
{
FatalErrorInFunction
<< "Problem: fast pointNeighbours routine did not find "
<< nbs.toc() << abort(FatalError);
}
}
}
// ************************************************************************* //

View File

@ -263,20 +263,35 @@ void Foam::patchDataWave<TransferType, TrackingData>::correct()
{
Map<label> nearestFace(2 * nWalls);
// Get distance and indices of nearest face
correctBoundaryFaceCells
(
patchIDs_,
distance_,
nearestFace
);
if (cellDistFuncs::useCombinedWallPatch)
{
// Correct across multiple patches
correctBoundaryCells
(
patchIDs_.sortedToc(),
true, // do point-connected cells as well
distance_,
nearestFace
);
}
else
{
// Get distance and indices of nearest face
correctBoundaryFaceCells
(
patchIDs_,
distance_,
nearestFace
);
correctBoundaryPointCells
(
patchIDs_,
distance_,
nearestFace
);
}
correctBoundaryPointCells
(
patchIDs_,
distance_,
nearestFace
);
// Transfer data from nearest face to cell
const List<TransferType>& faceInfo = waveInfo.allFaceInfo();

View File

@ -204,19 +204,34 @@ void Foam::patchWave::correct()
{
Map<label> nearestFace(2*nPatch);
correctBoundaryFaceCells
(
patchIDs_,
distance_,
nearestFace
);
if (cellDistFuncs::useCombinedWallPatch)
{
// Correct across multiple patches
correctBoundaryCells
(
patchIDs_.sortedToc(),
true, // do point-connected cells as well
distance_,
nearestFace
);
}
else
{
// Backwards compatible
correctBoundaryFaceCells
(
patchIDs_,
distance_,
nearestFace
);
correctBoundaryPointCells
(
patchIDs_,
distance_,
nearestFace
);
correctBoundaryPointCells
(
patchIDs_,
distance_,
nearestFace
);
}
}
}