zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: reduce boundary face queries in streamFunction

Browse Source 
BUG: streamFunction used uninitialized values for symmetry patches

- related to 8a8b5db977 changes (#3144)

ENH: improve robustness of surface field flattening

- vtk::surfaceFieldWriter
This commit is contained in:
Mark Olesen
2024-05-23 09:55:06 +02:00
parent 3f1d181b42
commit 6a80d4de40
5 changed files with 229 additions and 220 deletions
Download Patch File Download Diff File Expand all files Collapse all files

36
src/conversion/ensight/output/ensightOutputVolFieldTemplates.C
View File

@ -49,7 +49,7 @@ bool Foam::ensightOutput::writeVolField
bool parallel = Pstream::parRun(); bool parallel = Pstream::parRun();
const fvMesh& mesh = vf.mesh(); const fvMesh& mesh = vf.mesh();
const polyBoundaryMesh& bmesh = mesh.boundaryMesh(); const polyBoundaryMesh& pbm = mesh.boundaryMesh();
const Map<ensightCells>& cellZoneParts = ensMesh.cellZoneParts(); const Map<ensightCells>& cellZoneParts = ensMesh.cellZoneParts();
const Map<ensightFaces>& faceZoneParts = ensMesh.faceZoneParts(); const Map<ensightFaces>& faceZoneParts = ensMesh.faceZoneParts();
@ -69,13 +69,13 @@ bool Foam::ensightOutput::writeVolField
{ {
const ensightFaces& part = boundaryParts[patchId]; const ensightFaces& part = boundaryParts[patchId];
if (patchId < 0 || patchId >= bmesh.size()) if (patchId < 0 || patchId >= pbm.size())
{ {
// Future handling of combined patches? // Future handling of combined patches?
continue; continue;
} }
const label patchStart = bmesh[patchId].start(); const label patchStart = pbm[patchId].start();
// Either use a flat boundary field for all patches, // Either use a flat boundary field for all patches,
// or patch-local face ids // or patch-local face ids
@ -107,34 +107,28 @@ bool Foam::ensightOutput::writeVolField
// Flat boundary field // Flat boundary field
// similar to volPointInterpolation::flatBoundaryField() // 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) forAll(vf.boundaryField(), patchi)
{ {
const polyPatch& pp = bm[patchi].patch(); const polyPatch& pp = pbm[patchi];
const auto& bf = vf.boundaryField()[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 // Use average value for processor faces
// own cell value = patchInternalField // own cell value = patchInternalField
// nei cell value = evaluated boundary values // nei cell value = evaluated boundary values
SubList<Type> slice = (0.5 * (pfld.patchInternalField() + pfld));
(
flat,
bf.size(),
pp.offset()
) = (0.5 * (bf.patchInternalField() + bf));
} }
else if (!isA<emptyFvPatch>(bm[patchi])) else if (!isA<emptyPolyPatch>(pp))
{ {
SubList<Type> slice = pfld;
(
flat,
bf.size(),
pp.offset()
) = bf;
} }
} }

33
src/conversion/vtk/output/foamVtkSurfaceFieldWriter.C
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com \\ / A nd | www.openfoam.com
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
Copyright (C) 2016-2021 OpenCFD Ltd. Copyright (C) 2016-2024 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -26,35 +26,40 @@ License
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#include "foamVtkSurfaceFieldWriter.H" #include "foamVtkSurfaceFieldWriter.H"
#include "emptyFvsPatchFields.H"
#include "fvsPatchFields.H" #include "fvsPatchFields.H"
#include "surfaceFields.H" #include "surfaceFields.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
Foam::List<Foam::vector> Foam::vtk::surfaceFieldWriter::flattenBoundary namespace Foam
(
const surfaceVectorField& field
) const
{ {
// Boundary field - flatten
List<vector> flat(mesh_.nBoundaryFaces(), Zero); // Flatten boundary field values into a contiguous list
template<class Type>
static List<Type> flattenBoundary
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& field
)
{
const polyBoundaryMesh& pbm = field.mesh().boundaryMesh();
List<Type> flat(pbm.nFaces(), Foam::zero{});
forAll(field.boundaryField(), patchi) forAll(field.boundaryField(), patchi)
{ {
const polyPatch& pp = mesh_.boundaryMesh()[patchi]; const polyPatch& pp = pbm[patchi];
const auto& pfld = field.boundaryField()[patchi]; const auto& pfld = field.boundaryField()[patchi];
if (!isA<emptyFvsPatchVectorField>(pfld)) // Note: restrict transcribing to actual size of the patch field
{ // - handles "empty" patch type etc.
SubList<vector>(flat, pp.size(), pp.offset()) = pfld; SubList<Type>(flat, pfld.size(), pp.offset()) = pfld;
}
} }
return flat; return flat;
} }
} // End namespace Foam
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

8
src/conversion/vtk/output/foamVtkSurfaceFieldWriter.H
View File

@ -82,11 +82,9 @@ class surfaceFieldWriter
label numberOfPoints_; label numberOfPoints_;
// Private Member Functions public:
//- Flatten boundary field values into a contiguous list
List<vector> flattenBoundary(const surfaceVectorField& field) const;
// Generated Methods
//- No copy construct //- No copy construct
surfaceFieldWriter(const surfaceFieldWriter&) = delete; surfaceFieldWriter(const surfaceFieldWriter&) = delete;
@ -95,8 +93,6 @@ class surfaceFieldWriter
void operator=(const surfaceFieldWriter&) = delete; void operator=(const surfaceFieldWriter&) = delete;
public:
// Constructors // Constructors
//- Construct from mesh (default format INLINE_BASE64) //- Construct from mesh (default format INLINE_BASE64)

38
src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C
View File

@ -224,42 +224,34 @@ void Foam::volPointInterpolation::interpolateDimensionedInternalField
template<class Type> 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 GeometricField<Type, fvPatchField, volMesh>& vf
) const ) const
{ {
const fvMesh& mesh = vf.mesh(); const polyBoundaryMesh& pbm = vf.mesh().boundaryMesh();
const fvBoundaryMesh& bm = mesh.boundary();
auto tboundaryVals = tmp<Field<Type>>::New(mesh.nBoundaryFaces()); auto tboundaryVals = tmp<Field<Type>>::New(pbm.nFaces(), Foam::zero{});
auto& boundaryVals = tboundaryVals.ref(); auto& values = tboundaryVals.ref();
forAll(vf.boundaryField(), patchi) 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 if
( (
!isA<emptyFvPatch>(bm[patchi]) !isA<emptyPolyPatch>(pp)
&& !vf.boundaryField()[patchi].coupled() && !pfld.coupled()
) )
{ {
SubList<Type> slice = pfld;
(
boundaryVals,
vf.boundaryField()[patchi].size(),
bFacei
) = vf.boundaryField()[patchi];
}
else
{
const polyPatch& pp = bm[patchi].patch();
forAll(pp, i)
{
boundaryVals[bFacei++] = Zero;
}
} }
} }

334
src/functionObjects/field/streamFunction/streamFunction.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
Copyright (C) 2016 OpenFOAM Foundation Copyright (C) 2016 OpenFOAM Foundation
Copyright (C) 2019-2023 OpenCFD Ltd. Copyright (C) 2019-2024 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -77,9 +77,10 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
pMesh, pMesh,
dimensionedScalar(phi.dimensions(), Zero) 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 nVisited = 0;
label nVisitedOld = 0; label nVisitedOld = 0;
@ -87,10 +88,10 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
const faceUList& faces = mesh_.faces(); const faceUList& faces = mesh_.faces();
const pointField& points = mesh_.points(); const pointField& points = mesh_.points();
label nInternalFaces = mesh_.nInternalFaces(); const label nInternalFaces = mesh_.nInternalFaces();
vectorField unitAreas(mesh_.faceAreas()); vectorField unitAreas(mesh_.faceAreas());
unitAreas /= mag(unitAreas); unitAreas.normalise();
const polyPatchList& patches = mesh_.boundaryMesh(); const polyPatchList& patches = mesh_.boundaryMesh();
@ -104,45 +105,58 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
{ {
found = false; found = false;
// Check boundary faces first
forAll(patches, patchi) forAll(patches, patchi)
{ {
const primitivePatch& bouFaces = patches[patchi]; const auto& pp = patches[patchi];
const auto& patchPhi = phi.boundaryField()[patchi];
if (!isType<emptyPolyPatch>(patches[patchi])) // Skip empty, symmetry patches etc
if
(
patchPhi.empty()
|| isType<emptyPolyPatch>(pp)
|| isType<symmetryPlanePolyPatch>(pp)
|| isType<symmetryPolyPatch>(pp)
|| isType<wedgePolyPatch>(pp)
)
{ {
forAll(bouFaces, facei) continue;
}
forAll(pp, facei)
{
const auto& f = pp[facei];
if (magSqr(patchPhi[facei]) < SMALL)
{ {
if (magSqr(phi.boundaryField()[patchi][facei]) < SMALL) // Zero flux face found
found = true;
for (const label pointi : f)
{ {
const labelList& zeroPoints = bouFaces[facei]; if (visitedPoint.test(pointi))
// Zero flux face found
found = true;
forAll(zeroPoints, pointi)
{ {
if (visitedPoint[zeroPoints[pointi]] == 1) found = false;
{
found = false;
break;
}
}
if (found)
{
Log << " Zero face: patch: " << patchi
<< " face: " << facei << endl;
forAll(zeroPoints, pointi)
{
streamFunction[zeroPoints[pointi]] = 0;
visitedPoint[zeroPoints[pointi]] = 1;
nVisited++;
}
break; break;
} }
} }
if (found)
{
Log << " Zero face: patch: " << patchi
<< " face: " << facei << endl;
for (const label pointi : f)
{
visitedPoint.set(pointi);
++nVisited;
streamFunction[pointi] = 0;
}
break;
}
} }
} }
@ -152,20 +166,17 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
if (!found) if (!found)
{ {
Log << " Zero flux boundary face not found. " Log << " Zero flux boundary face not found. "
<< "Using cell as a reference." << "Using cell as a reference." << endl;
<< endl;
const cellList& c = mesh_.cells(); for (const cell& c : mesh_.cells())
forAll(c, ci)
{ {
labelList zeroPoints = c[ci].labels(mesh_.faces()); labelList zeroPoints = c.labels(mesh_.faces());
bool found = true; bool found = true;
forAll(zeroPoints, pointi) for (const label pointi : zeroPoints)
{ {
if (visitedPoint[zeroPoints[pointi]] == 1) if (visitedPoint.test(pointi))
{ {
found = false; found = false;
break; break;
@ -174,11 +185,12 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
if (found) if (found)
{ {
forAll(zeroPoints, pointi) for (const label pointi : zeroPoints)
{ {
streamFunction[zeroPoints[pointi]] = 0.0; visitedPoint.set(pointi);
visitedPoint[zeroPoints[pointi]] = 1; ++nVisited;
nVisited++;
streamFunction[pointi] = 0;
} }
break; break;
@ -201,132 +213,136 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
{ {
finished = true; finished = true;
for (label facei = nInternalFaces; facei<faces.size(); facei++) scalar currentStreamValue(0);
point currentStreamPoint(Zero);
// Boundary faces first
forAll(patches, patchi)
{ {
const labelList& curBPoints = faces[facei]; const auto& pp = patches[patchi];
bool bPointFound = false; const auto& patchPhi = phi.boundaryField()[patchi];
scalar currentBStream = 0.0; // Skip empty, symmetry patches etc
vector currentBStreamPoint(0, 0, 0); if
(
forAll(curBPoints, pointi) patchPhi.empty()
|| isType<emptyPolyPatch>(pp)
|| isType<symmetryPlanePolyPatch>(pp)
|| isType<symmetryPolyPatch>(pp)
|| isType<wedgePolyPatch>(pp)
)
{ {
// Check if the point has been visited continue;
if (visitedPoint[curBPoints[pointi]] == 1)
{
// The point has been visited
currentBStream = streamFunction[curBPoints[pointi]];
currentBStreamPoint = points[curBPoints[pointi]];
bPointFound = true;
break;
}
} }
if (bPointFound) forAll(pp, facei)
{ {
// Sort out other points on the face const auto& f = pp[facei];
forAll(curBPoints, pointi)
// Check if the point has been visited
bool pointFound = false;
for (const label pointi : f)
{ {
// Check if the point has been visited if (visitedPoint.test(pointi))
if (visitedPoint[curBPoints[pointi]] == 0)
{ {
label patchNo = // The point has been visited
mesh_.boundaryMesh().whichPatch(facei); currentStreamValue = streamFunction[pointi];
currentStreamPoint = points[pointi];
if pointFound = true;
break;
}
}
if (!pointFound)
{
finished = false;
continue;
}
// Sort out other points on the face
for (const label pointi : f)
{
// If the point has not yet been visited
if (!visitedPoint.test(pointi))
{
vector edgeHat =
( (
!isType<emptyPolyPatch>(patches[patchNo]) points[pointi] - currentStreamPoint
&& !isType<symmetryPlanePolyPatch> );
(patches[patchNo]) edgeHat.replace(slabDir, 0);
&& !isType<symmetryPolyPatch>(patches[patchNo]) edgeHat.normalise();
&& !isType<wedgePolyPatch>(patches[patchNo])
) const vector& nHat = unitAreas[facei];
if (edgeHat.y() > VSMALL)
{ {
label faceNo = visitedPoint.set(pointi);
mesh_.boundaryMesh()[patchNo] ++nVisited;
.whichFace(facei);
vector edgeHat = streamFunction[pointi] =
points[curBPoints[pointi]] (
- currentBStreamPoint; currentStreamValue
edgeHat.replace(slabDir, 0); + patchPhi[facei]*sign(nHat.x())
edgeHat.normalise(); );
}
else if (edgeHat.y() < -VSMALL)
{
visitedPoint.set(pointi);
++nVisited;
vector nHat = unitAreas[facei]; streamFunction[pointi] =
(
if (edgeHat.y() > VSMALL) currentStreamValue
- patchPhi[facei]*sign(nHat.x())
);
}
else
{
if (edgeHat.x() > VSMALL)
{ {
visitedPoint[curBPoints[pointi]] = 1; visitedPoint.set(pointi);
nVisited++; ++nVisited;
streamFunction[curBPoints[pointi]] = streamFunction[pointi] =
currentBStream (
+ phi.boundaryField()[patchNo][faceNo] currentStreamValue
*sign(nHat.x()); + patchPhi[facei]*sign(nHat.y())
);
} }
else if (edgeHat.y() < -VSMALL) else if (edgeHat.x() < -VSMALL)
{ {
visitedPoint[curBPoints[pointi]] = 1; visitedPoint.set(pointi);
nVisited++; ++nVisited;
streamFunction[curBPoints[pointi]] = streamFunction[pointi] =
currentBStream (
- phi.boundaryField()[patchNo][faceNo] currentStreamValue
*sign(nHat.x()); - patchPhi[facei]*sign(nHat.y())
} );
else
{
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());
}
} }
} }
} }
} }
} }
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 auto& f = faces[facei];
const labelList& curPoints = faces[facei];
bool pointFound = false; bool pointFound = false;
scalar currentStream = 0.0; for (const label pointi : f)
point currentStreamPoint(0, 0, 0);
forAll(curPoints, pointi)
{ {
// Check if the point has been visited // Check if the point has been visited
if (visitedPoint[curPoints[pointi]] == 1) if (visitedPoint.test(pointi))
{ {
// The point has been visited currentStreamValue = streamFunction[pointi];
currentStream = streamFunction[curPoints[pointi]]; currentStreamPoint = points[pointi];
currentStreamPoint = points[curPoints[pointi]];
pointFound = true; pointFound = true;
break; break;
@ -336,36 +352,42 @@ Foam::tmp<Foam::pointScalarField> Foam::functionObjects::streamFunction::calc
if (pointFound) if (pointFound)
{ {
// Sort out other points on the face // Sort out other points on the face
forAll(curPoints, pointi) for (const label pointi : f)
{ {
// Check if the point has been visited // If the point has not yet been visited
if (visitedPoint[curPoints[pointi]] == 0) if (!visitedPoint.test(pointi))
{ {
vector edgeHat = vector edgeHat =
points[curPoints[pointi]] - currentStreamPoint; (
points[pointi] - currentStreamPoint
);
edgeHat.replace(slabDir, 0); edgeHat.replace(slabDir, 0);
edgeHat.normalise(); edgeHat.normalise();
vector nHat = unitAreas[facei]; const vector& nHat = unitAreas[facei];
if (edgeHat.y() > VSMALL) if (edgeHat.y() > VSMALL)
{ {
visitedPoint[curPoints[pointi]] = 1; visitedPoint.set(pointi);
nVisited++; ++nVisited;
streamFunction[curPoints[pointi]] = streamFunction[pointi] =
currentStream (
+ phi[facei]*sign(nHat.x()); currentStreamValue
+ phi[facei]*sign(nHat.x())
);
} }
else if (edgeHat.y() < -VSMALL) else if (edgeHat.y() < -VSMALL)
{ {
visitedPoint[curPoints[pointi]] = 1; visitedPoint.set(pointi);
nVisited++; ++nVisited;
streamFunction[curPoints[pointi]] = streamFunction[pointi] =
currentStream (
- phi[facei]*sign(nHat.x()); 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(); const scalar thickness = vector(slabNormal) & mesh_.bounds().span();
streamFunction /= thickness; streamFunction /= thickness;
streamFunction.boundaryFieldRef() = 0.0; streamFunction.boundaryFieldRef() = Zero;
return tstreamFunction; return tstreamFunction;
} }