diff --git a/src/conversion/ensight/output/ensightOutputVolFieldTemplates.C b/src/conversion/ensight/output/ensightOutputVolFieldTemplates.C index b940f11d23..0a2cb59146 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& cellZoneParts = ensMesh.cellZoneParts(); const Map& 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 flat(mesh.nBoundaryFaces(), Zero); + Field 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(bm[patchi])) + // Note: restrict transcribing to actual size of the patch field + // - handles "empty" patch type etc. + + SubList slice(flat, pfld.size(), pp.offset()); + + if (isA(pp)) { // Use average value for processor faces // own cell value = patchInternalField // nei cell value = evaluated boundary values - SubList - ( - flat, - bf.size(), - pp.offset() - ) = (0.5 * (bf.patchInternalField() + bf)); + slice = (0.5 * (pfld.patchInternalField() + pfld)); } - else if (!isA(bm[patchi])) + else if (!isA(pp)) { - SubList - ( - 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 8bb1a1ee49..c6494bd880 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::vtk::surfaceFieldWriter::flattenBoundary -( - const surfaceVectorField& field -) const +namespace Foam { - // Boundary field - flatten - List flat(mesh_.nBoundaryFaces(), Zero); +// Flatten boundary field values into a contiguous list +template +static List flattenBoundary +( + const GeometricField& field +) +{ + const polyBoundaryMesh& pbm = field.mesh().boundaryMesh(); + + List 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(pfld)) - { - SubList(flat, pp.size(), pp.offset()) = pfld; - } + // Note: restrict transcribing to actual size of the patch field + // - handles "empty" patch type etc. + SubList(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 1dc5ca1c6d..420fb0fca6 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 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 fac2318f22..3ef984cf4d 100644 --- a/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C +++ b/src/finiteVolume/interpolation/volPointInterpolation/volPointInterpolate.C @@ -224,42 +224,34 @@ void Foam::volPointInterpolation::interpolateDimensionedInternalField template -Foam::tmp> Foam::volPointInterpolation::flatBoundaryField +Foam::tmp> +Foam::volPointInterpolation::flatBoundaryField ( const GeometricField& vf ) const { - const fvMesh& mesh = vf.mesh(); - const fvBoundaryMesh& bm = mesh.boundary(); + const polyBoundaryMesh& pbm = vf.mesh().boundaryMesh(); - auto tboundaryVals = tmp>::New(mesh.nBoundaryFaces()); - auto& boundaryVals = tboundaryVals.ref(); + auto tboundaryVals = tmp>::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 slice(values, pfld.size(), pp.offset()); if ( - !isA(bm[patchi]) - && !vf.boundaryField()[patchi].coupled() + !isA(pp) + && !pfld.coupled() ) { - SubList - ( - 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 741efa28d4..9380b1591d 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::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::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,45 +105,58 @@ Foam::tmp 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]; - if (!isType(patches[patchi])) + // Skip empty, symmetry patches etc + if + ( + patchPhi.empty() + || isType(pp) + || isType(pp) + || isType(pp) + || isType(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]; - - // Zero flux face found - found = true; - - forAll(zeroPoints, pointi) + if (visitedPoint.test(pointi)) { - if (visitedPoint[zeroPoints[pointi]] == 1) - { - 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++; - } - + found = false; 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::functionObjects::streamFunction::calc if (!found) { Log << " Zero flux boundary face not found. " - << "Using cell as a reference." - << endl; + << "Using cell as a reference." << endl; - const cellList& c = mesh_.cells(); - - 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::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::functionObjects::streamFunction::calc { finished = true; - for (label facei = nInternalFaces; facei(pp) + || isType(pp) + || isType(pp) + || isType(pp) + ) { - // Check if the point has been visited - if (visitedPoint[curBPoints[pointi]] == 1) - { - // The point has been visited - currentBStream = streamFunction[curBPoints[pointi]]; - currentBStreamPoint = points[curBPoints[pointi]]; - - bPointFound = true; - - break; - } + continue; } - if (bPointFound) + forAll(pp, facei) { - // Sort out other points on the face - forAll(curBPoints, pointi) + const auto& f = pp[facei]; + + // Check if the point has been visited + bool pointFound = false; + + for (const label pointi : f) { - // Check if the point has been visited - if (visitedPoint[curBPoints[pointi]] == 0) + if (visitedPoint.test(pointi)) { - label patchNo = - mesh_.boundaryMesh().whichPatch(facei); + // The point has been visited + 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(patches[patchNo]) - && !isType - (patches[patchNo]) - && !isType(patches[patchNo]) - && !isType(patches[patchNo]) - ) + points[pointi] - currentStreamPoint + ); + edgeHat.replace(slabDir, 0); + edgeHat.normalise(); + + const vector& nHat = unitAreas[facei]; + + if (edgeHat.y() > VSMALL) { - label faceNo = - mesh_.boundaryMesh()[patchNo] - .whichFace(facei); + visitedPoint.set(pointi); + ++nVisited; - vector edgeHat = - points[curBPoints[pointi]] - - currentBStreamPoint; - edgeHat.replace(slabDir, 0); - edgeHat.normalise(); + streamFunction[pointi] = + ( + currentStreamValue + + patchPhi[facei]*sign(nHat.x()) + ); + } + else if (edgeHat.y() < -VSMALL) + { + visitedPoint.set(pointi); + ++nVisited; - vector nHat = unitAreas[facei]; - - if (edgeHat.y() > VSMALL) + streamFunction[pointi] = + ( + currentStreamValue + - patchPhi[facei]*sign(nHat.x()) + ); + } + else + { + if (edgeHat.x() > VSMALL) { - visitedPoint[curBPoints[pointi]] = 1; - nVisited++; + visitedPoint.set(pointi); + ++nVisited; - streamFunction[curBPoints[pointi]] = - currentBStream - + phi.boundaryField()[patchNo][faceNo] - *sign(nHat.x()); + streamFunction[pointi] = + ( + currentStreamValue + + patchPhi[facei]*sign(nHat.y()) + ); } - else if (edgeHat.y() < -VSMALL) + else if (edgeHat.x() < -VSMALL) { - visitedPoint[curBPoints[pointi]] = 1; - nVisited++; + visitedPoint.set(pointi); + ++nVisited; - streamFunction[curBPoints[pointi]] = - currentBStream - - phi.boundaryField()[patchNo][faceNo] - *sign(nHat.x()); - } - 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()); - } + streamFunction[pointi] = + ( + currentStreamValue + - patchPhi[facei]*sign(nHat.y()) + ); } } } } } - else - { - finished = false; - } } - for (label facei=0; facei 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++; + visitedPoint.set(pointi); + ++nVisited; - streamFunction[curPoints[pointi]] = - currentStream - + phi[facei]*sign(nHat.x()); + streamFunction[pointi] = + ( + currentStreamValue + + phi[facei]*sign(nHat.x()) + ); } else if (edgeHat.y() < -VSMALL) { - visitedPoint[curPoints[pointi]] = 1; - nVisited++; + visitedPoint.set(pointi); + ++nVisited; - streamFunction[curPoints[pointi]] = - currentStream - - phi[facei]*sign(nHat.x()); + streamFunction[pointi] = + ( + currentStreamValue + - phi[facei]*sign(nHat.x()) + ); } } } @@ -397,7 +419,7 @@ Foam::tmp Foam::functionObjects::streamFunction::calc const scalar thickness = vector(slabNormal) & mesh_.bounds().span(); streamFunction /= thickness; - streamFunction.boundaryFieldRef() = 0.0; + streamFunction.boundaryFieldRef() = Zero; return tstreamFunction; }