From 6a80d4de400da755f44fd95189af6aa45c5551dd Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Thu, 23 May 2024 09:55:06 +0200
Subject: [PATCH] ENH: reduce boundary face queries in streamFunction
BUG: streamFunction used uninitialized values for symmetry patches
- related to 8a8b5db977c3 changes (#3144)
ENH: improve robustness of surface field flattening
- vtk::surfaceFieldWriter
---
.../output/ensightOutputVolFieldTemplates.C | 36 +-
.../vtk/output/foamVtkSurfaceFieldWriter.C | 31 +-
.../vtk/output/foamVtkSurfaceFieldWriter.H | 8 +-
.../volPointInterpolate.C | 38 +-
.../field/streamFunction/streamFunction.C | 332 ++++++++++--------
5 files changed, 227 insertions(+), 218 deletions(-)
diff --git a/src/conversion/ensight/output/ensightOutputVolFieldTemplates.C b/src/conversion/ensight/output/ensightOutputVolFieldTemplates.C
index b940f11d23b..0a2cb59146e 100644
--- a/src/conversion/ensight/output/ensightOutputVolFieldTemplates.C
+++ b/src/conversion/ensight/output/ensightOutputVolFieldTemplates.C
@@ -49,7 +49,7 @@ bool Foam::ensightOutput::writeVolField
bool parallel = Pstream::parRun();
const fvMesh& mesh = vf.mesh();
- const polyBoundaryMesh& bmesh = mesh.boundaryMesh();
+ const polyBoundaryMesh& pbm = mesh.boundaryMesh();
const Map<ensightCells>& cellZoneParts = ensMesh.cellZoneParts();
const Map<ensightFaces>& faceZoneParts = ensMesh.faceZoneParts();
@@ -69,13 +69,13 @@ bool Foam::ensightOutput::writeVolField
{
const ensightFaces& part = boundaryParts[patchId];
- if (patchId < 0 || patchId >= bmesh.size())
+ if (patchId < 0 || patchId >= pbm.size())
{
// Future handling of combined patches?
continue;
}
- const label patchStart = bmesh[patchId].start();
+ const label patchStart = pbm[patchId].start();
// Either use a flat boundary field for all patches,
// or patch-local face ids
@@ -107,34 +107,28 @@ bool Foam::ensightOutput::writeVolField
// Flat boundary field
// similar to volPointInterpolation::flatBoundaryField()
- Field<Type> flat(mesh.nBoundaryFaces(), Zero);
+ Field<Type> flat(pbm.nFaces(), Foam::zero{});
- const fvBoundaryMesh& bm = mesh.boundary();
forAll(vf.boundaryField(), patchi)
{
- const polyPatch& pp = bm[patchi].patch();
- const auto& bf = vf.boundaryField()[patchi];
+ const polyPatch& pp = pbm[patchi];
+ const auto& pfld = vf.boundaryField()[patchi];
- if (isA<processorFvPatch>(bm[patchi]))
+ // Note: restrict transcribing to actual size of the patch field
+ // - handles "empty" patch type etc.
+
+ SubList<Type> slice(flat, pfld.size(), pp.offset());
+
+ if (isA<processorPolyPatch>(pp))
{
// Use average value for processor faces
// own cell value = patchInternalField
// nei cell value = evaluated boundary values
- SubList<Type>
- (
- flat,
- bf.size(),
- pp.offset()
- ) = (0.5 * (bf.patchInternalField() + bf));
+ slice = (0.5 * (pfld.patchInternalField() + pfld));
}
- else if (!isA<emptyFvPatch>(bm[patchi]))
+ else if (!isA<emptyPolyPatch>(pp))
{
- SubList<Type>
- (
- flat,
- bf.size(),
- pp.offset()
- ) = bf;
+ slice = pfld;
}
}
diff --git a/src/conversion/vtk/output/foamVtkSurfaceFieldWriter.C b/src/conversion/vtk/output/foamVtkSurfaceFieldWriter.C
index 8bb1a1ee49e..c6494bd880d 100644
--- a/src/conversion/vtk/output/foamVtkSurfaceFieldWriter.C
+++ b/src/conversion/vtk/output/foamVtkSurfaceFieldWriter.C
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2016-2021 OpenCFD Ltd.
+ Copyright (C) 2016-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -26,35 +26,40 @@ License
\*---------------------------------------------------------------------------*/
#include "foamVtkSurfaceFieldWriter.H"
-#include "emptyFvsPatchFields.H"
#include "fvsPatchFields.H"
#include "surfaceFields.H"
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
-Foam::List<Foam::vector> Foam::vtk::surfaceFieldWriter::flattenBoundary
+namespace Foam
+{
+
+// Flatten boundary field values into a contiguous list
+template<class Type>
+static List<Type> flattenBoundary
(
- const surfaceVectorField& field
-) const
+ const GeometricField<Type, fvsPatchField, surfaceMesh>& field
+)
{
- // Boundary field - flatten
+ const polyBoundaryMesh& pbm = field.mesh().boundaryMesh();
- List<vector> flat(mesh_.nBoundaryFaces(), Zero);
+ List<Type> flat(pbm.nFaces(), Foam::zero{});
forAll(field.boundaryField(), patchi)
{
- const polyPatch& pp = mesh_.boundaryMesh()[patchi];
+ const polyPatch& pp = pbm[patchi];
const auto& pfld = field.boundaryField()[patchi];
- if (!isA<emptyFvsPatchVectorField>(pfld))
- {
- SubList<vector>(flat, pp.size(), pp.offset()) = pfld;
- }
+ // Note: restrict transcribing to actual size of the patch field
+ // - handles "empty" patch type etc.
+ SubList<Type>(flat, pfld.size(), pp.offset()) = pfld;
}
return flat;
}
+} // End namespace Foam
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
diff --git a/src/conversion/vtk/output/foamVtkSurfaceFieldWriter.H b/src/conversion/vtk/output/foamVtkSurfaceFieldWriter.H
index 1dc5ca1c6d7..420fb0fca6d 100644
--- a/src/conversion/vtk/output/foamVtkSurfaceFieldWriter.H
+++ b/src/conversion/vtk/output/foamVtkSurfaceFieldWriter.H
@@ -82,11 +82,9 @@ class surfaceFieldWriter
label numberOfPoints_;
- // Private Member Functions
-
- //- Flatten boundary field values into a contiguous list
- List<vector> flattenBoundary(const surfaceVectorField& field) const;
+public:
+ // Generated Methods
//- No copy construct
surfaceFieldWriter(const surfaceFieldWriter&) = delete;
@@ -95,8 +93,6 @@ class surfaceFieldWriter
void operator=(const surfaceFieldWriter&) = delete;
-public:
-
// Constructors
//- Construct from mesh (default format INLINE_BASE64)
diff --git a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
index fac2318f22d..3ef984cf4d7 100644
--- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
+++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
@@ -224,42 +224,34 @@ void Foam::volPointInterpolation::interpolateDimensionedInternalField
template<class Type>
-Foam::tmp<Foam::Field<Type>> Foam::volPointInterpolation::flatBoundaryField
+Foam::tmp<Foam::Field<Type>>
+Foam::volPointInterpolation::flatBoundaryField
(
const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
- const fvMesh& mesh = vf.mesh();
- const fvBoundaryMesh& bm = mesh.boundary();
+ const polyBoundaryMesh& pbm = vf.mesh().boundaryMesh();
- auto tboundaryVals = tmp<Field<Type>>::New(mesh.nBoundaryFaces());
- auto& boundaryVals = tboundaryVals.ref();
+ auto tboundaryVals = tmp<Field<Type>>::New(pbm.nFaces(), Foam::zero{});
+ auto& values = tboundaryVals.ref();
forAll(vf.boundaryField(), patchi)
{
- label bFacei = bm[patchi].patch().start() - mesh.nInternalFaces();
+ const auto& pp = pbm[patchi];
+ const auto& pfld = vf.boundaryField()[patchi];
+
+ // Note: restrict transcribing to actual size of the patch field
+ // - handles "empty" patch type etc.
+
+ SubList<Type> slice(values, pfld.size(), pp.offset());
if
(
- !isA<emptyFvPatch>(bm[patchi])
- && !vf.boundaryField()[patchi].coupled()
+ !isA<emptyPolyPatch>(pp)
+ && !pfld.coupled()
)
{
- SubList<Type>
- (
- boundaryVals,
- vf.boundaryField()[patchi].size(),
- bFacei
- ) = vf.boundaryField()[patchi];
- }
- else
- {
- const polyPatch& pp = bm[patchi].patch();
-
- forAll(pp, i)
- {
- boundaryVals[bFacei++] = Zero;
- }
+ slice = pfld;
}
}
diff --git a/src/functionObjects/field/streamFunction/streamFunction.C b/src/functionObjects/field/streamFunction/streamFunction.C
index 741efa28d4a..9380b1591d0 100644
--- a/src/functionObjects/field/streamFunction/streamFunction.C
+++ b/src/functionObjects/field/streamFunction/streamFunction.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation
- Copyright (C) 2019-2023 OpenCFD Ltd.
+ Copyright (C) 2019-2024 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -77,9 +77,10 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
pMesh,
dimensionedScalar(phi.dimensions(), Zero)
);
- pointScalarField& streamFunction = tstreamFunction.ref();
+ auto& streamFunction = tstreamFunction.ref();
- labelList visitedPoint(mesh_.nPoints(), Zero);
+
+ bitSet visitedPoint(mesh_.nPoints());
label nVisited = 0;
label nVisitedOld = 0;
@@ -87,10 +88,10 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
const faceUList& faces = mesh_.faces();
const pointField& points = mesh_.points();
- label nInternalFaces = mesh_.nInternalFaces();
+ const label nInternalFaces = mesh_.nInternalFaces();
vectorField unitAreas(mesh_.faceAreas());
- unitAreas /= mag(unitAreas);
+ unitAreas.normalise();
const polyPatchList& patches = mesh_.boundaryMesh();
@@ -104,44 +105,57 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
{
found = false;
+ // Check boundary faces first
forAll(patches, patchi)
{
- const primitivePatch& bouFaces = patches[patchi];
+ const auto& pp = patches[patchi];
+ const auto& patchPhi = phi.boundaryField()[patchi];
+
+ // Skip empty, symmetry patches etc
+ if
+ (
+ patchPhi.empty()
+ || isType<emptyPolyPatch>(pp)
+ || isType<symmetryPlanePolyPatch>(pp)
+ || isType<symmetryPolyPatch>(pp)
+ || isType<wedgePolyPatch>(pp)
+ )
+ {
+ continue;
+ }
- if (!isType<emptyPolyPatch>(patches[patchi]))
+ forAll(pp, facei)
{
- forAll(bouFaces, facei)
- {
- if (magSqr(phi.boundaryField()[patchi][facei]) < SMALL)
- {
- const labelList& zeroPoints = bouFaces[facei];
+ const auto& f = pp[facei];
- // Zero flux face found
- found = true;
+ if (magSqr(patchPhi[facei]) < SMALL)
+ {
+ // Zero flux face found
+ found = true;
- forAll(zeroPoints, pointi)
+ for (const label pointi : f)
+ {
+ if (visitedPoint.test(pointi))
{
- if (visitedPoint[zeroPoints[pointi]] == 1)
- {
- found = false;
- break;
- }
+ found = false;
+ break;
}
+ }
- if (found)
- {
- Log << " Zero face: patch: " << patchi
- << " face: " << facei << endl;
+ if (found)
+ {
+ Log << " Zero face: patch: " << patchi
+ << " face: " << facei << endl;
- forAll(zeroPoints, pointi)
- {
- streamFunction[zeroPoints[pointi]] = 0;
- visitedPoint[zeroPoints[pointi]] = 1;
- nVisited++;
- }
+ for (const label pointi : f)
+ {
+ visitedPoint.set(pointi);
+ ++nVisited;
- break;
+ streamFunction[pointi] = 0;
}
+
+ break;
}
}
}
@@ -152,20 +166,17 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
if (!found)
{
Log << " Zero flux boundary face not found. "
- << "Using cell as a reference."
- << endl;
-
- const cellList& c = mesh_.cells();
+ << "Using cell as a reference." << endl;
- forAll(c, ci)
+ for (const cell& c : mesh_.cells())
{
- labelList zeroPoints = c[ci].labels(mesh_.faces());
+ labelList zeroPoints = c.labels(mesh_.faces());
bool found = true;
- forAll(zeroPoints, pointi)
+ for (const label pointi : zeroPoints)
{
- if (visitedPoint[zeroPoints[pointi]] == 1)
+ if (visitedPoint.test(pointi))
{
found = false;
break;
@@ -174,11 +185,12 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
if (found)
{
- forAll(zeroPoints, pointi)
+ for (const label pointi : zeroPoints)
{
- streamFunction[zeroPoints[pointi]] = 0.0;
- visitedPoint[zeroPoints[pointi]] = 1;
- nVisited++;
+ visitedPoint.set(pointi);
+ ++nVisited;
+
+ streamFunction[pointi] = 0;
}
break;
@@ -201,132 +213,136 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
{
finished = true;
- for (label facei = nInternalFaces; facei<faces.size(); facei++)
- {
- const labelList& curBPoints = faces[facei];
- bool bPointFound = false;
+ scalar currentStreamValue(0);
+ point currentStreamPoint(Zero);
- scalar currentBStream = 0.0;
- vector currentBStreamPoint(0, 0, 0);
+ // Boundary faces first
+ forAll(patches, patchi)
+ {
+ const auto& pp = patches[patchi];
+ const auto& patchPhi = phi.boundaryField()[patchi];
+
+ // Skip empty, symmetry patches etc
+ if
+ (
+ patchPhi.empty()
+ || isType<emptyPolyPatch>(pp)
+ || isType<symmetryPlanePolyPatch>(pp)
+ || isType<symmetryPolyPatch>(pp)
+ || isType<wedgePolyPatch>(pp)
+ )
+ {
+ continue;
+ }
- forAll(curBPoints, pointi)
+ forAll(pp, facei)
{
+ const auto& f = pp[facei];
+
// Check if the point has been visited
- if (visitedPoint[curBPoints[pointi]] == 1)
+ bool pointFound = false;
+
+ for (const label pointi : f)
{
- // The point has been visited
- currentBStream = streamFunction[curBPoints[pointi]];
- currentBStreamPoint = points[curBPoints[pointi]];
+ if (visitedPoint.test(pointi))
+ {
+ // The point has been visited
+ currentStreamValue = streamFunction[pointi];
+ currentStreamPoint = points[pointi];
- bPointFound = true;
+ pointFound = true;
+ break;
+ }
+ }
- break;
+ if (!pointFound)
+ {
+ finished = false;
+ continue;
}
- }
- if (bPointFound)
- {
+
// Sort out other points on the face
- forAll(curBPoints, pointi)
+ for (const label pointi : f)
{
- // Check if the point has been visited
- if (visitedPoint[curBPoints[pointi]] == 0)
+ // If the point has not yet been visited
+ if (!visitedPoint.test(pointi))
{
- label patchNo =
- mesh_.boundaryMesh().whichPatch(facei);
-
- if
+ vector edgeHat =
(
- !isType<emptyPolyPatch>(patches[patchNo])
- && !isType<symmetryPlanePolyPatch>
- (patches[patchNo])
- && !isType<symmetryPolyPatch>(patches[patchNo])
- && !isType<wedgePolyPatch>(patches[patchNo])
- )
- {
- label faceNo =
- mesh_.boundaryMesh()[patchNo]
- .whichFace(facei);
-
- vector edgeHat =
- points[curBPoints[pointi]]
- - currentBStreamPoint;
- edgeHat.replace(slabDir, 0);
- edgeHat.normalise();
+ points[pointi] - currentStreamPoint
+ );
+ edgeHat.replace(slabDir, 0);
+ edgeHat.normalise();
- vector nHat = unitAreas[facei];
+ const vector& nHat = unitAreas[facei];
- if (edgeHat.y() > VSMALL)
- {
- visitedPoint[curBPoints[pointi]] = 1;
- nVisited++;
-
- streamFunction[curBPoints[pointi]] =
- currentBStream
- + phi.boundaryField()[patchNo][faceNo]
- *sign(nHat.x());
- }
- else if (edgeHat.y() < -VSMALL)
+ if (edgeHat.y() > VSMALL)
+ {
+ visitedPoint.set(pointi);
+ ++nVisited;
+
+ streamFunction[pointi] =
+ (
+ currentStreamValue
+ + patchPhi[facei]*sign(nHat.x())
+ );
+ }
+ else if (edgeHat.y() < -VSMALL)
+ {
+ visitedPoint.set(pointi);
+ ++nVisited;
+
+ streamFunction[pointi] =
+ (
+ currentStreamValue
+ - patchPhi[facei]*sign(nHat.x())
+ );
+ }
+ else
+ {
+ if (edgeHat.x() > VSMALL)
{
- visitedPoint[curBPoints[pointi]] = 1;
- nVisited++;
-
- streamFunction[curBPoints[pointi]] =
- currentBStream
- - phi.boundaryField()[patchNo][faceNo]
- *sign(nHat.x());
+ visitedPoint.set(pointi);
+ ++nVisited;
+
+ streamFunction[pointi] =
+ (
+ currentStreamValue
+ + patchPhi[facei]*sign(nHat.y())
+ );
}
- else
+ else if (edgeHat.x() < -VSMALL)
{
- if (edgeHat.x() > VSMALL)
- {
- visitedPoint[curBPoints[pointi]] = 1;
- nVisited++;
-
- streamFunction[curBPoints[pointi]] =
- currentBStream
- + phi.boundaryField()[patchNo][faceNo]
- *sign(nHat.y());
- }
- else if (edgeHat.x() < -VSMALL)
- {
- visitedPoint[curBPoints[pointi]] = 1;
- nVisited++;
-
- streamFunction[curBPoints[pointi]] =
- currentBStream
- - phi.boundaryField()[patchNo][faceNo]
- *sign(nHat.y());
- }
+ visitedPoint.set(pointi);
+ ++nVisited;
+
+ streamFunction[pointi] =
+ (
+ currentStreamValue
+ - patchPhi[facei]*sign(nHat.y())
+ );
}
}
}
}
}
- else
- {
- finished = false;
- }
}
- for (label facei=0; facei<nInternalFaces; facei++)
+ // Internal faces next
+ for (label facei = 0; facei < nInternalFaces; ++facei)
{
- // Get the list of point labels for the face
- const labelList& curPoints = faces[facei];
+ const auto& f = faces[facei];
bool pointFound = false;
- scalar currentStream = 0.0;
- point currentStreamPoint(0, 0, 0);
-
- forAll(curPoints, pointi)
+ for (const label pointi : f)
{
// Check if the point has been visited
- if (visitedPoint[curPoints[pointi]] == 1)
+ if (visitedPoint.test(pointi))
{
- // The point has been visited
- currentStream = streamFunction[curPoints[pointi]];
- currentStreamPoint = points[curPoints[pointi]];
+ currentStreamValue = streamFunction[pointi];
+ currentStreamPoint = points[pointi];
pointFound = true;
break;
@@ -336,36 +352,42 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
if (pointFound)
{
// Sort out other points on the face
- forAll(curPoints, pointi)
+ for (const label pointi : f)
{
- // Check if the point has been visited
- if (visitedPoint[curPoints[pointi]] == 0)
+ // If the point has not yet been visited
+ if (!visitedPoint.test(pointi))
{
vector edgeHat =
- points[curPoints[pointi]] - currentStreamPoint;
+ (
+ points[pointi] - currentStreamPoint
+ );
edgeHat.replace(slabDir, 0);
edgeHat.normalise();
- vector nHat = unitAreas[facei];
+ const vector& nHat = unitAreas[facei];
if (edgeHat.y() > VSMALL)
{
- visitedPoint[curPoints[pointi]] = 1;
- nVisited++;
-
- streamFunction[curPoints[pointi]] =
- currentStream
- + phi[facei]*sign(nHat.x());
+ visitedPoint.set(pointi);
+ ++nVisited;
+
+ streamFunction[pointi] =
+ (
+ currentStreamValue
+ + phi[facei]*sign(nHat.x())
+ );
}
else if (edgeHat.y() < -VSMALL)
{
- visitedPoint[curPoints[pointi]] = 1;
- nVisited++;
-
- streamFunction[curPoints[pointi]] =
- currentStream
- - phi[facei]*sign(nHat.x());
+ visitedPoint.set(pointi);
+ ++nVisited;
+
+ streamFunction[pointi] =
+ (
+ currentStreamValue
+ - phi[facei]*sign(nHat.x())
+ );
}
}
}
@@ -397,7 +419,7 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
const scalar thickness = vector(slabNormal) & mesh_.bounds().span();
streamFunction /= thickness;
- streamFunction.boundaryFieldRef() = 0.0;
+ streamFunction.boundaryFieldRef() = Zero;
return tstreamFunction;
}
--
GitLab