From 60b31fc8e29f0e5850eded92573551ffb842ee29 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Mon, 4 Apr 2022 18:48:17 +0200 Subject: [PATCH] ENH: add primitiveMeshTools support for face lists - allows reuse by finiteArea, for example. - simplify edge looping with face thisLabel/nextLabel method ENH: additional storage checks for mesh weights (faMesh + fvMesh) - allow finite-area field decomposition without edge weights. STYLE: use tmp New in various places. Simpler updateGeom check --- .../polyMesh/polyMeshCheck/polyMeshTools.C | 45 +++-- .../polyMesh/polyMeshCheck/polyMeshTools.H | 17 +- .../PatchTools/PatchToolsNormals.C | 4 +- .../meshes/primitiveMesh/primitiveMesh.C | 25 ++- .../meshes/primitiveMesh/primitiveMesh.H | 8 +- .../primitiveMeshCellCentresAndVols.C | 9 +- .../primitiveMeshCheck/primitiveMeshTools.C | 163 +++++++++++------- .../primitiveMeshCheck/primitiveMeshTools.H | 75 ++++++-- .../meshes/primitiveMesh/primitiveMeshClear.C | 13 +- .../primitiveMeshFaceCentresAndAreas.C | 9 +- .../meshes/primitiveMesh/primitiveMeshI.H | 12 +- .../polyMeshGeometry/polyMeshGeometry.C | 56 +++--- .../polyTopoChange/combineFaces.C | 15 +- .../edgeInterpolation/edgeInterpolation.C | 6 +- .../edgeInterpolation/edgeInterpolation.H | 52 +++--- .../surfaceInterpolation.H | 15 +- .../primitiveMeshGeometry.C | 7 +- .../decompose/faDecompose/faFieldDecomposer.C | 12 +- 18 files changed, 313 insertions(+), 230 deletions(-) diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C index a90525c3bb..6f03bfb165 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2012-2016 OpenFOAM Foundation - Copyright (C) 2021 OpenCFD Ltd. + Copyright (C) 2021-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -31,7 +31,7 @@ License #include "pyramidPointFaceRef.H" #include "primitiveMeshTools.H" -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // Foam::tmp Foam::polyMeshTools::faceOrthogonality ( @@ -44,8 +44,8 @@ Foam::tmp Foam::polyMeshTools::faceOrthogonality const labelList& nei = mesh.faceNeighbour(); const polyBoundaryMesh& pbm = mesh.boundaryMesh(); - tmp tortho(new scalarField(mesh.nFaces(), 1.0)); - scalarField& ortho = tortho.ref(); + auto tortho = tmp::New(mesh.nFaces(), scalar(1)); + auto& ortho = tortho.ref(); // Internal faces forAll(nei, facei) @@ -64,9 +64,8 @@ Foam::tmp Foam::polyMeshTools::faceOrthogonality pointField neighbourCc; syncTools::swapBoundaryCellPositions(mesh, cc, neighbourCc); - forAll(pbm, patchi) + for (const polyPatch& pp : pbm) { - const polyPatch& pp = pbm[patchi]; if (pp.coupled()) { forAll(pp, i) @@ -99,16 +98,17 @@ Foam::tmp Foam::polyMeshTools::faceSkewness { const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); + const faceList& faces = mesh.faces(); const polyBoundaryMesh& pbm = mesh.boundaryMesh(); - tmp tskew(new scalarField(mesh.nFaces())); - scalarField& skew = tskew.ref(); + auto tskew = tmp::New(mesh.nFaces()); + auto& skew = tskew.ref(); forAll(nei, facei) { skew[facei] = primitiveMeshTools::faceSkewness ( - mesh, + faces, p, fCtrs, fAreas, @@ -126,9 +126,8 @@ Foam::tmp Foam::polyMeshTools::faceSkewness pointField neighbourCc; syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neighbourCc); - forAll(pbm, patchi) + for (const polyPatch& pp : pbm) { - const polyPatch& pp = pbm[patchi]; if (pp.coupled()) { forAll(pp, i) @@ -138,7 +137,7 @@ Foam::tmp Foam::polyMeshTools::faceSkewness skew[facei] = primitiveMeshTools::faceSkewness ( - mesh, + faces, p, fCtrs, fAreas, @@ -157,7 +156,7 @@ Foam::tmp Foam::polyMeshTools::faceSkewness skew[facei] = primitiveMeshTools::boundaryFaceSkewness ( - mesh, + faces, p, fCtrs, fAreas, @@ -185,8 +184,8 @@ Foam::tmp Foam::polyMeshTools::faceWeights const labelList& nei = mesh.faceNeighbour(); const polyBoundaryMesh& pbm = mesh.boundaryMesh(); - tmp tweight(new scalarField(mesh.nFaces(), 1.0)); - scalarField& weight = tweight.ref(); + auto tweight = tmp::New(mesh.nFaces(), scalar(1)); + auto& weight = tweight.ref(); // Internal faces forAll(nei, facei) @@ -206,9 +205,8 @@ Foam::tmp Foam::polyMeshTools::faceWeights pointField neiCc; syncTools::swapBoundaryCellPositions(mesh, cellCtrs, neiCc); - forAll(pbm, patchi) + for (const polyPatch& pp : pbm) { - const polyPatch& pp = pbm[patchi]; if (pp.coupled()) { forAll(pp, i) @@ -241,8 +239,8 @@ Foam::tmp Foam::polyMeshTools::volRatio const labelList& nei = mesh.faceNeighbour(); const polyBoundaryMesh& pbm = mesh.boundaryMesh(); - tmp tratio(new scalarField(mesh.nFaces(), 1.0)); - scalarField& ratio = tratio.ref(); + auto tratio = tmp::New(mesh.nFaces(), scalar(1)); + auto& ratio = tratio.ref(); // Internal faces forAll(nei, facei) @@ -259,9 +257,8 @@ Foam::tmp Foam::polyMeshTools::volRatio scalarField neiVol; syncTools::swapBoundaryCellList(mesh, vol, neiVol); - forAll(pbm, patchi) + for (const polyPatch& pp : pbm) { - const polyPatch& pp = pbm[patchi]; if (pp.coupled()) { forAll(pp, i) @@ -305,10 +302,8 @@ Foam::polyMesh::readUpdateState Foam::polyMeshTools::combine { return state1; } - else - { - return state0; - } + + return state0; } diff --git a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H index 8e599d9e75..a7da2ec716 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H +++ b/src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.H @@ -24,7 +24,7 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . -Namespace +Class Foam::polyMeshTools Description @@ -36,8 +36,8 @@ SourceFiles \*---------------------------------------------------------------------------*/ -#ifndef polyMeshTools_H -#define polyMeshTools_H +#ifndef Foam_polyMeshTools_H +#define Foam_polyMeshTools_H #include "polyMesh.H" #include "primitiveMeshTools.H" @@ -48,18 +48,17 @@ namespace Foam { /*---------------------------------------------------------------------------*\ - Namespace polyMeshTools Declaration + Class polyMeshTools Declaration \*---------------------------------------------------------------------------*/ class polyMeshTools : public primitiveMeshTools { - public: - //- Generate orthogonality field. (1 for fully orthogonal, < 1 for - // non-orthogonal) + //- Generate orthogonality field. + //- (1 for fully orthogonal, < 1 for non-orthogonal) static tmp faceOrthogonality ( const polyMesh& mesh, @@ -93,8 +92,8 @@ public: const scalarField& vol ); - //- Combine readUpdateState. topo change trumps geom-only - // change etc. + //- Combine readUpdateState. + //- topo change trumps geom-only change etc. static polyMesh::readUpdateState combine ( const polyMesh::readUpdateState& state0, diff --git a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C index f606cbc50d..f99a940bf9 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C +++ b/src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.C @@ -156,8 +156,8 @@ Foam::PatchTools::pointNormals // 1. Start off with local normals (note:without calculating pointNormals // to avoid them being stored) - tmp textrudeN(new pointField(p.nPoints(), Zero)); - pointField& extrudeN = textrudeN.ref(); + auto textrudeN = tmp::New(p.nPoints(), Zero); + auto& extrudeN = textrudeN.ref(); { const faceList& localFaces = p.localFaces(); const vectorField& faceNormals = p.faceNormals(); diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C index 1b468aa414..997a8113bd 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.C @@ -333,14 +333,14 @@ Foam::tmp Foam::primitiveMesh::movePoints } // Create swept volumes - const faceList& f = faces(); + const faceList& fcs = faces(); - tmp tsweptVols(new scalarField(f.size())); - scalarField& sweptVols = tsweptVols.ref(); + auto tsweptVols = tmp::New(fcs.size()); + auto& sweptVols = tsweptVols.ref(); - forAll(f, facei) + forAll(fcs, facei) { - sweptVols[facei] = f[facei].sweptVol(oldPoints, newPoints); + sweptVols[facei] = fcs[facei].sweptVol(oldPoints, newPoints); } // Force recalculation of all geometric data with new points @@ -364,20 +364,15 @@ const Foam::cellShapeList& Foam::primitiveMesh::cellShapes() const void Foam::primitiveMesh::updateGeom() { - if (!faceCentresPtr_) + if (!faceCentresPtr_ || !faceAreasPtr_) { + // These are always calculated in tandem, but only once calcFaceCentresAndAreas(); } - if (!faceAreasPtr_) - { - calcFaceCentresAndAreas(); - } - if (!cellCentresPtr_) - { - calcCellCentresAndVols(); - } - if (!cellVolumesPtr_) + + if (!cellCentresPtr_ || !cellVolumesPtr_) { + // These are always calculated in tandem, but only once calcCellCentresAndVols(); } } diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H index f807041e49..5bcda5e24d 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMesh.H @@ -50,8 +50,8 @@ SourceFiles \*---------------------------------------------------------------------------*/ -#ifndef primitiveMesh_H -#define primitiveMesh_H +#ifndef Foam_primitiveMesh_H +#define Foam_primitiveMesh_H #include "DynamicList.H" #include "edgeList.H" @@ -188,7 +188,7 @@ class primitiveMesh void operator=(const primitiveMesh&) = delete; - // Topological calculations + // Topological Calculations //- Calculate cell shapes void calcCellShapes() const; @@ -809,8 +809,8 @@ public: inline bool hasPointPoints() const noexcept; inline bool hasCellPoints() const noexcept; inline bool hasCellCentres() const noexcept; - inline bool hasFaceCentres() const noexcept; inline bool hasCellVolumes() const noexcept; + inline bool hasFaceCentres() const noexcept; inline bool hasFaceAreas() const noexcept; // On-the-fly addressing calculation. These functions return either diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C index f8a1e32340..71eb30e7b9 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCellCentresAndVols.C @@ -39,16 +39,15 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const if (debug) { Pout<< "primitiveMesh::calcCellCentresAndVols() : " - << "Calculating cell centres and cell volumes" + << "Calculating cell centres and volumes" << endl; } - // It is an error to attempt to recalculate cellCentres - // if the pointer is already set + // These are always calculated in tandem, but only once. if (cellCentresPtr_ || cellVolumesPtr_) { FatalErrorInFunction - << "Cell centres or cell volumes already calculated" + << "Cell centres or volumes already calculated" << abort(FatalError); } @@ -73,7 +72,7 @@ void Foam::primitiveMesh::calcCellCentresAndVols() const if (debug) { Pout<< "primitiveMesh::calcCellCentresAndVols() : " - << "Finished calculating cell centres and cell volumes" + << "Finished calculating cell centres and volumes" << endl; } } diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C index 9694b933a8..c2e5d1cdd8 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2012-2016 OpenFOAM Foundation - Copyright (C) 2017-2021 OpenCFD Ltd. + Copyright (C) 2017-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -32,21 +32,23 @@ License #include "pyramidPointFaceRef.H" #include "PrecisionAdaptor.H" -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // void Foam::primitiveMeshTools::makeFaceCentresAndAreas ( - const primitiveMesh& mesh, + const UList& fcs, const pointField& p, vectorField& fCtrs, vectorField& fAreas ) { - const faceList& fs = mesh.faces(); + // Safety first - ensure properly sized + fCtrs.resize_nocopy(fcs.size()); + fAreas.resize_nocopy(fcs.size()); - forAll(fs, facei) + forAll(fcs, facei) { - const labelList& f = fs[facei]; + const face& f = fcs[facei]; const label nPoints = f.size(); // If the face is a triangle, do a direct calculation for efficiency @@ -59,26 +61,25 @@ void Foam::primitiveMeshTools::makeFaceCentresAndAreas else { solveVector sumN = Zero; - solveScalar sumA = 0.0; + solveScalar sumA = Zero; solveVector sumAc = Zero; solveVector fCentre = p[f[0]]; - for (label pi = 1; pi < nPoints; pi++) + for (label pi = 1; pi < nPoints; ++pi) { fCentre += solveVector(p[f[pi]]); } - fCentre /= nPoints; - for (label pi = 0; pi < nPoints; pi++) + for (label pi = 0; pi < nPoints; ++pi) { - const label nextPi(pi == nPoints-1 ? 0 : pi+1); - const solveVector nextPoint(p[f[nextPi]]); - const solveVector thisPoint(p[f[pi]]); + const solveVector thisPoint(p[f.thisLabel(pi)]); + const solveVector nextPoint(p[f.nextLabel(pi)]); solveVector c = thisPoint + nextPoint + fCentre; solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint); solveScalar a = mag(n); + sumN += n; sumA += a; sumAc += a*c; @@ -101,6 +102,18 @@ void Foam::primitiveMeshTools::makeFaceCentresAndAreas } +void Foam::primitiveMeshTools::makeFaceCentresAndAreas +( + const primitiveMesh& mesh, + const pointField& p, + vectorField& fCtrs, + vectorField& fAreas +) +{ + makeFaceCentresAndAreas(mesh.faces(), p, fCtrs, fAreas); +} + + void Foam::primitiveMeshTools::makeCellCentresAndVols ( const primitiveMesh& mesh, @@ -117,7 +130,7 @@ void Foam::primitiveMeshTools::makeCellCentresAndVols // Clear the fields for accumulation cellCtrs = Zero; - cellVols = 0.0; + cellVols = Zero; const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); @@ -201,7 +214,7 @@ void Foam::primitiveMeshTools::makeCellCentresAndVols Foam::scalar Foam::primitiveMeshTools::faceSkewness ( - const primitiveMesh& mesh, + const UList& fcs, const pointField& p, const vectorField& fCtrs, const vectorField& fAreas, @@ -211,6 +224,8 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness const point& neiCc ) { + const face& f = fcs[facei]; + vector Cpf = fCtrs[facei] - ownCc; vector d = neiCc - ownCc; @@ -224,7 +239,6 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness // from the face centre to the edge of the face in the direction // of the skewness scalar fd = 0.2*mag(d) + ROOTVSMALL; - const face& f = mesh.faces()[facei]; forAll(f, pi) { fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei]))); @@ -235,6 +249,60 @@ Foam::scalar Foam::primitiveMeshTools::faceSkewness } +Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness +( + const UList& fcs, + const pointField& p, + const vectorField& fCtrs, + const vectorField& fAreas, + + const label facei, + const point& ownCc +) +{ + const face& f = fcs[facei]; + + vector Cpf = fCtrs[facei] - ownCc; + + vector normal = normalised(fAreas[facei]); + vector d = normal*(normal & Cpf); + + // Skewness vector + vector sv = + Cpf + - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d; + vector svHat = sv/(mag(sv) + ROOTVSMALL); + + // Normalisation distance calculated as the approximate distance + // from the face centre to the edge of the face in the direction + // of the skewness + scalar fd = 0.4*mag(d) + ROOTVSMALL; + forAll(f, pi) + { + fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei]))); + } + + // Normalised skewness + return mag(sv)/fd; +} + + +Foam::scalar Foam::primitiveMeshTools::faceSkewness +( + const primitiveMesh& mesh, + const pointField& p, + const vectorField& fCtrs, + const vectorField& fAreas, + + const label facei, + const point& ownCc, + const point& neiCc +) +{ + return faceSkewness(mesh.faces(), p, fCtrs, fAreas, facei, ownCc, neiCc); +} + + Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness ( const primitiveMesh& mesh, @@ -246,30 +314,7 @@ Foam::scalar Foam::primitiveMeshTools::boundaryFaceSkewness const point& ownCc ) { - vector Cpf = fCtrs[facei] - ownCc; - - vector normal = normalised(fAreas[facei]); - vector d = normal*(normal & Cpf); - - - // Skewness vector - vector sv = - Cpf - - ((fAreas[facei] & Cpf)/((fAreas[facei] & d) + ROOTVSMALL))*d; - vector svHat = sv/(mag(sv) + ROOTVSMALL); - - // Normalisation distance calculated as the approximate distance - // from the face centre to the edge of the face in the direction - // of the skewness - scalar fd = 0.4*mag(d) + ROOTVSMALL; - const face& f = mesh.faces()[facei]; - forAll(f, pi) - { - fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[facei]))); - } - - // Normalised skewness - return mag(sv)/fd; + return boundaryFaceSkewness(mesh.faces(), p, fCtrs, fAreas, facei, ownCc); } @@ -298,8 +343,8 @@ Foam::tmp Foam::primitiveMeshTools::faceOrthogonality const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); - tmp tortho(new scalarField(mesh.nInternalFaces())); - scalarField& ortho = tortho.ref(); + auto tortho = tmp::New(mesh.nInternalFaces()); + auto& ortho = tortho.ref(); // Internal faces forAll(nei, facei) @@ -327,15 +372,17 @@ Foam::tmp Foam::primitiveMeshTools::faceSkewness { const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); + const faceList& faces = mesh.faces(); - tmp tskew(new scalarField(mesh.nFaces())); - scalarField& skew = tskew.ref(); + auto tskew = tmp::New(mesh.nFaces()); + auto& skew = tskew.ref(); - forAll(nei, facei) + // Internal faces + for (label facei = 0; facei < mesh.nInternalFaces(); ++facei) { skew[facei] = faceSkewness ( - mesh, + faces, p, fCtrs, fAreas, @@ -350,11 +397,11 @@ Foam::tmp Foam::primitiveMeshTools::faceSkewness // Boundary faces: consider them to have only skewness error. // (i.e. treat as if mirror cell on other side) - for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); facei++) + for (label facei = mesh.nInternalFaces(); facei < mesh.nFaces(); ++facei) { skew[facei] = boundaryFaceSkewness ( - mesh, + faces, p, fCtrs, fAreas, @@ -513,15 +560,14 @@ Foam::tmp Foam::primitiveMeshTools::faceConcavity vectorField faceNormals(faceAreas); faceNormals /= mag(faceNormals) + ROOTVSMALL; - tmp tfaceAngles(new scalarField(mesh.nFaces())); - scalarField& faceAngles = tfaceAngles.ref(); - + auto tfaceAngles = tmp::New(mesh.nFaces()); + auto&& faceAngles = tfaceAngles.ref(); forAll(fcs, facei) { const face& f = fcs[facei]; - // Get edge from f[0] to f[size-1]; + // Normalized vector from f[size-1] to f[0]; vector ePrev(p[f.first()] - p[f.last()]); scalar magEPrev = mag(ePrev); ePrev /= magEPrev + ROOTVSMALL; @@ -530,11 +576,8 @@ Foam::tmp Foam::primitiveMeshTools::faceConcavity forAll(f, fp0) { - // Get vertex after fp - label fp1 = f.fcIndex(fp0); - // Normalized vector between two consecutive points - vector e10(p[f[fp1]] - p[f[fp0]]); + vector e10(p[f.nextLabel(fp0)] - p[f.thisLabel(fp0)]); scalar magE10 = mag(e10); e10 /= magE10 + ROOTVSMALL; @@ -584,8 +627,8 @@ Foam::tmp Foam::primitiveMeshTools::faceFlatness // primitiveMeshFaceCentresAndAreas.C) scalarField magAreas(mag(faceAreas)); - tmp tfaceFlatness(new scalarField(mesh.nFaces(), 1.0)); - scalarField& faceFlatness = tfaceFlatness.ref(); + auto tfaceFlatness = tmp::New(mesh.nFaces(), scalar(1)); + auto& faceFlatness = tfaceFlatness.ref(); forAll(fcs, facei) { @@ -641,8 +684,8 @@ Foam::tmp Foam::primitiveMeshTools::cellDeterminant } } - tmp tcellDeterminant(new scalarField(mesh.nCells())); - scalarField& cellDeterminant = tcellDeterminant.ref(); + auto tcellDeterminant = tmp::New(mesh.nCells()); + auto& cellDeterminant = tcellDeterminant.ref(); const cellList& c = mesh.cells(); diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H index 3bfef13593..b5001a3322 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2012-2016 OpenFOAM Foundation - Copyright (C) 2017 OpenCFD Ltd. + Copyright (C) 2017-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -24,7 +24,7 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . -Namespace +Class Foam::primitiveMeshTools Description @@ -35,8 +35,8 @@ SourceFiles \*---------------------------------------------------------------------------*/ -#ifndef primitiveMeshTools_H -#define primitiveMeshTools_H +#ifndef Foam_primitiveMeshTools_H +#define Foam_primitiveMeshTools_H #include "bitSet.H" #include "primitiveFields.H" @@ -47,24 +47,36 @@ SourceFiles namespace Foam { -// Forward declarations +// Forward Declarations +class face; class primitiveMesh; /*---------------------------------------------------------------------------*\ - Namespace primitiveMeshTools Declaration + Class primitiveMeshTools Declaration \*---------------------------------------------------------------------------*/ class primitiveMeshTools { public: - //- Calculate face centres and areas + //- Calculate face centres and areas for specified faces. + // Adjusts the lengths of centres and area normals if required. static void makeFaceCentresAndAreas ( - const primitiveMesh& mesh, - const pointField& p, - vectorField& fCtrs, - vectorField& fAreas + const UList& faces, //!< The faces to consider + const pointField& p, //!< Support points for faces + vectorField& fCtrs, //!< [out] The face centres + vectorField& fAreas //!< [out] The face area normals + ); + + //- Calculate face centres and areas for all mesh faces. + // Adjusts the lengths of centres and area normals if required. + static void makeFaceCentresAndAreas + ( + const primitiveMesh& mesh, //!< The mesh faces + const pointField& p, //!< Support points for faces + vectorField& fCtrs, //!< [out] The face centres + vectorField& fAreas //!< [out] The face area normals ); //- Calculate cell centres and volumes from face properties @@ -117,7 +129,7 @@ public: ); //- Generate face concavity field. Returns per face the (sin of the) - // most concave angle between two consecutive edges + //- most concave angle between two consecutive edges static tmp faceConcavity ( const scalar maxSin, @@ -127,8 +139,8 @@ public: ); //- Generate face flatness field. Compares the individual triangles' - // normals against the face average normal. Between 0 (fully warped) - // and 1 (fully flat) + //- normals against the face average normal. Between 0 (fully warped) + //- and 1 (fully flat) static tmp faceFlatness ( const primitiveMesh& mesh, @@ -138,7 +150,7 @@ public: ); //- Generate edge alignment field. Is per face the minimum aligned edge - // (does not use edge addressing) + //- (does not use edge addressing) static tmp edgeAlignment ( const primitiveMesh& mesh, @@ -161,8 +173,8 @@ public: //- Skewness of single face static scalar faceSkewness ( - const primitiveMesh& mesh, - const pointField& p, + const UList& faces, //!< The faces to consider + const pointField& p, //!< Support points for faces const vectorField& fCtrs, const vectorField& fAreas, @@ -174,8 +186,33 @@ public: //- Skewness of single boundary face static scalar boundaryFaceSkewness ( - const primitiveMesh& mesh, - const pointField& p, + const UList& faces, //!< The faces to consider + const pointField& p, //!< Support points for faces + const vectorField& fCtrs, + const vectorField& fAreas, + + const label facei, + const point& ownCc + ); + + //- Skewness of single face + static scalar faceSkewness + ( + const primitiveMesh& mesh, //!< The mesh faces + const pointField& p, //!< Support points for faces + const vectorField& fCtrs, + const vectorField& fAreas, + + const label facei, + const point& ownCc, + const point& neiCc + ); + + //- Skewness of single boundary face + static scalar boundaryFaceSkewness + ( + const primitiveMesh& mesh, //!< The mesh faces + const pointField& p, //!< Support points for faces const vectorField& fCtrs, const vectorField& fAreas, diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C index bf19ed1ffe..13494c0cae 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshClear.C @@ -106,21 +106,20 @@ void Foam::primitiveMesh::printAllocated() const Pout<< " Cell-centres" << endl; } - if (faceCentresPtr_) - { - Pout<< " Face-centres" << endl; - } - if (cellVolumesPtr_) { Pout<< " Cell-volumes" << endl; } + if (faceCentresPtr_) + { + Pout<< " Face-centres" << endl; + } + if (faceAreasPtr_) { Pout<< " Face-areas" << endl; } - } @@ -134,8 +133,8 @@ void Foam::primitiveMesh::clearGeom() } deleteDemandDrivenData(cellCentresPtr_); - deleteDemandDrivenData(faceCentresPtr_); deleteDemandDrivenData(cellVolumesPtr_); + deleteDemandDrivenData(faceCentresPtr_); deleteDemandDrivenData(faceAreasPtr_); } diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C index dda3e3634c..f3577d04ec 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshFaceCentresAndAreas.C @@ -42,16 +42,15 @@ void Foam::primitiveMesh::calcFaceCentresAndAreas() const if (debug) { Pout<< "primitiveMesh::calcFaceCentresAndAreas() : " - << "Calculating face centres and face areas" + << "Calculating face centres and areas" << endl; } - // It is an error to attempt to recalculate faceCentres - // if the pointer is already set + // These are always calculated in tandem, but only once. if (faceCentresPtr_ || faceAreasPtr_) { FatalErrorInFunction - << "Face centres or face areas already calculated" + << "Face centres and areas already calculated" << abort(FatalError); } @@ -66,7 +65,7 @@ void Foam::primitiveMesh::calcFaceCentresAndAreas() const if (debug) { Pout<< "primitiveMesh::calcFaceCentresAndAreas() : " - << "Finished calculating face centres and face areas" + << "Finished calculating face centres and areas" << endl; } } diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H index e3e70b2394..1a877dfbf3 100644 --- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H +++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshI.H @@ -192,18 +192,18 @@ inline bool Foam::primitiveMesh::hasCellCentres() const noexcept } -inline bool Foam::primitiveMesh::hasFaceCentres() const noexcept -{ - return bool(faceCentresPtr_); -} - - inline bool Foam::primitiveMesh::hasCellVolumes() const noexcept { return bool(cellVolumesPtr_); } +inline bool Foam::primitiveMesh::hasFaceCentres() const noexcept +{ + return bool(faceCentresPtr_); +} + + inline bool Foam::primitiveMesh::hasFaceAreas() const noexcept { return bool(faceAreasPtr_); diff --git a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C index 06fac222cc..59a2cce429 100644 --- a/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C +++ b/src/dynamicMesh/motionSmoother/polyMeshGeometry/polyMeshGeometry.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2017 OpenFOAM Foundation - Copyright (C) 2015-2019 OpenCFD Ltd. + Copyright (C) 2015-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -50,12 +50,12 @@ void Foam::polyMeshGeometry::updateFaceCentresAndAreas const labelList& changedFaces ) { - const faceList& fs = mesh_.faces(); + const faceList& fcs = mesh_.faces(); for (const label facei : changedFaces) { - const labelList& f = fs[facei]; - label nPoints = f.size(); + const face& f = fcs[facei]; + const label nPoints = f.size(); // If the face is a triangle, do a direct calculation for efficiency // and to avoid round-off error-related problems @@ -66,33 +66,41 @@ void Foam::polyMeshGeometry::updateFaceCentresAndAreas } else { - vector sumN = Zero; - scalar sumA = Zero; - vector sumAc = Zero; + solveVector sumN = Zero; + solveScalar sumA = Zero; + solveVector sumAc = Zero; - point fCentre = p[f[0]]; + solveVector fCentre = p[f[0]]; for (label pi = 1; pi < nPoints; ++pi) { - fCentre += p[f[pi]]; + fCentre += solveVector(p[f[pi]]); } - fCentre /= nPoints; for (label pi = 0; pi < nPoints; ++pi) { - const point& nextPoint = p[f[(pi + 1) % nPoints]]; + const solveVector thisPoint(p[f.thisLabel(pi)]); + const solveVector nextPoint(p[f.nextLabel(pi)]); - vector c = p[f[pi]] + nextPoint + fCentre; - vector n = (nextPoint - p[f[pi]])^(fCentre - p[f[pi]]); - scalar a = mag(n); + solveVector c = thisPoint + nextPoint + fCentre; + solveVector n = (nextPoint - thisPoint)^(fCentre - thisPoint); + solveScalar a = mag(n); sumN += n; sumA += a; sumAc += a*c; } - faceCentres_[facei] = (1.0/3.0)*sumAc/(sumA + VSMALL); - faceAreas_[facei] = 0.5*sumN; + if (sumA < ROOTVSMALL) + { + faceCentres_[facei] = fCentre; + faceAreas_[facei] = Zero; + } + else + { + faceCentres_[facei] = (1.0/3.0)*sumAc/sumA; + faceAreas_[facei] = 0.5*sumN; + } } } } @@ -946,6 +954,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness const labelList& own = mesh.faceOwner(); const labelList& nei = mesh.faceNeighbour(); + const faceList& faces = mesh.faces(); const polyBoundaryMesh& patches = mesh.boundaryMesh(); // Calculate coupled cell centre @@ -962,7 +971,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness { scalar skewness = primitiveMeshTools::faceSkewness ( - mesh, + faces, points, faceCentres, faceAreas, @@ -996,7 +1005,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness { scalar skewness = primitiveMeshTools::faceSkewness ( - mesh, + faces, points, faceCentres, faceAreas, @@ -1030,7 +1039,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness { scalar skewness = primitiveMeshTools::boundaryFaceSkewness ( - mesh, + faces, points, faceCentres, faceAreas, @@ -1072,7 +1081,7 @@ bool Foam::polyMeshGeometry::checkFaceSkewness scalar skewness = primitiveMeshTools::faceSkewness ( - mesh, + faces, points, faceCentres, faceAreas, @@ -1449,18 +1458,15 @@ bool Foam::polyMeshGeometry::checkFaceAngles const vector faceNormal = normalised(faceAreas[facei]); - // Get edge from f[0] to f[size-1]; + // Normalized vector from f[size-1] to f[0]; vector ePrev(p[f.first()] - p[f.last()]); scalar magEPrev = mag(ePrev); ePrev /= magEPrev + VSMALL; forAll(f, fp0) { - // Get vertex after fp - label fp1 = f.fcIndex(fp0); - // Normalized vector between two consecutive points - vector e10(p[f[fp1]] - p[f[fp0]]); + vector e10(p[f.nextLabel(fp0)] - p[f.thisLabel(fp0)]); scalar magE10 = mag(e10); e10 /= magE10 + VSMALL; diff --git a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C index 02517d8f49..21e5944423 100644 --- a/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C +++ b/src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2018-2021 OpenCFD Ltd. + Copyright (C) 2018-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -61,18 +61,15 @@ bool Foam::combineFaces::convexFace // Get outwards pointing normal of f, only the sign matters. const vector areaNorm = f.areaNormal(points); - // Get edge from f[0] to f[size-1]; + // Normalized vector from f[size-1] to f[0]; vector ePrev(points[f.first()] - points[f.last()]); scalar magEPrev = mag(ePrev); ePrev /= magEPrev + VSMALL; forAll(f, fp0) { - // Get vertex after fp - label fp1 = f.fcIndex(fp0); - // Normalized vector between two consecutive points - vector e10(points[f[fp1]] - points[f[fp0]]); + vector e10(points[f.nextLabel(fp0)] - points[f.thisLabel(fp0)]); scalar magE10 = mag(e10); e10 /= magE10 + VSMALL; @@ -138,12 +135,10 @@ void Foam::combineFaces::regioniseFaces { const polyBoundaryMesh& patches = mesh_.boundaryMesh(); - forAll(cEdges, i) + for (const label edgei : cEdges) { - const label edgeI = cEdges[i]; - label f0, f1; - meshTools::getEdgeFaces(mesh_, celli, edgeI, f0, f1); + meshTools::getEdgeFaces(mesh_, celli, edgei, f0, f1); const vector& a0 = mesh_.faceAreas()[f0]; const vector& a1 = mesh_.faceAreas()[f1]; diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C index 0ed885bab5..370e95ff6d 100644 --- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C +++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.C @@ -59,10 +59,10 @@ Foam::edgeInterpolation::edgeInterpolation(const faMesh& fam) lPN_(nullptr), weightingFactors_(nullptr), differenceFactors_(nullptr), - orthogonal_(false), correctionVectors_(nullptr), - skew_(true), - skewCorrectionVectors_(nullptr) + skewCorrectionVectors_(nullptr), + orthogonal_(false), + skew_(true) {} diff --git a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H index 995256d9bb..b740b0a07f 100644 --- a/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H +++ b/src/finiteArea/interpolation/edgeInterpolation/edgeInterpolation.H @@ -38,8 +38,8 @@ SourceFiles \*---------------------------------------------------------------------------*/ -#ifndef edgeInterpolation_H -#define edgeInterpolation_H +#ifndef Foam_edgeInterpolation_H +#define Foam_edgeInterpolation_H #include "tmp.H" #include "scalar.H" @@ -52,6 +52,7 @@ SourceFiles namespace Foam { +// Forward Declarations class polyMesh; /*---------------------------------------------------------------------------*\ @@ -60,7 +61,7 @@ class polyMesh; class edgeInterpolation { - // Private data + // Private Data // Reference to faMesh const faMesh& faMesh_; @@ -77,20 +78,20 @@ class edgeInterpolation //- Face-gradient difference factors mutable edgeScalarField* differenceFactors_; - //- Is mesh orthogonal - mutable bool orthogonal_; - //- Non-orthogonality correction vectors mutable edgeVectorField* correctionVectors_; - //- Is mesh skew - mutable bool skew_; - //- Skew correction vectors mutable edgeVectorField* skewCorrectionVectors_; + //- Is mesh orthogonal + mutable bool orthogonal_; - // Private member functions + //- Is mesh skew + mutable bool skew_; + + + // Private Member Functions //- Construct geodesic distance between P and N void makeLPN() const; @@ -107,15 +108,12 @@ class edgeInterpolation //- Construct skewness correction vectors void makeSkewCorrectionVectors() const; -// //- Construct Least-squares gradient vectors -// void makeLeastSquareVectors() const; - protected: - // Protected member functions + // Protected Member Functions - // Storage management + // Storage Management //- Clear all geometry and addressing void clearOut(); @@ -137,10 +135,10 @@ public: ~edgeInterpolation(); - // Member functions + // Member Functions //- Return mesh reference - const faMesh& mesh() const + const faMesh& mesh() const noexcept { return faMesh_; } @@ -154,21 +152,29 @@ public: //- Return reference to difference factors array const edgeScalarField& deltaCoeffs() const; - //- Return whether mesh is orthogonal or not - bool orthogonal() const; - //- Return reference to non-orthogonality correction vectors array const edgeVectorField& correctionVectors() const; - //- Return whether mesh is skew or not - bool skew() const; - //- Return reference to skew vectors array const edgeVectorField& skewCorrectionVectors() const; + //- Return whether mesh is orthogonal or not + bool orthogonal() const; + + //- Return whether mesh is skew or not + bool skew() const; + + + // Mesh Motion //- Do what is necessary if the mesh has moved bool movePoints() const; + + + // Storage Management + + //- True if weights exist + bool hasWeights() const noexcept { return bool(weightingFactors_); } }; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H index 5f83411a19..03ab721be7 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H +++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.H @@ -35,8 +35,8 @@ SourceFiles \*---------------------------------------------------------------------------*/ -#ifndef surfaceInterpolation_H -#define surfaceInterpolation_H +#ifndef Foam_surfaceInterpolation_H +#define Foam_surfaceInterpolation_H #include "tmp.H" #include "scalar.H" @@ -49,6 +49,7 @@ SourceFiles namespace Foam { +// Forward Declarations class fvMesh; class fvGeometryScheme; @@ -58,7 +59,7 @@ class fvGeometryScheme; class surfaceInterpolation { - // Private data + // Private Data // Reference to fvMesh const fvMesh& mesh_; @@ -108,7 +109,7 @@ public: virtual ~surfaceInterpolation(); - // Member functions + // Member Functions //- Return reference to geometry calculation scheme virtual const fvGeometryScheme& geometry() const; @@ -134,6 +135,12 @@ public: //- Update all geometric data virtual void updateGeom(); + + + // Storage Management + + //- Has weights + bool hasWeights() const noexcept { return bool(weights_); } }; diff --git a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C index 9b063d22e6..5628472da8 100644 --- a/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C +++ b/src/meshTools/primitiveMeshGeometry/primitiveMeshGeometry.C @@ -763,18 +763,15 @@ bool Foam::primitiveMeshGeometry::checkFaceAngles const vector faceNormal = normalised(faceAreas[facei]); - // Get edge from f[0] to f[size-1]; + // Normalized vector from f[size-1] to f[0]; vector ePrev(p[f.first()] - p[f.last()]); scalar magEPrev = mag(ePrev); ePrev /= magEPrev + VSMALL; forAll(f, fp0) { - // Get vertex after fp - label fp1 = f.fcIndex(fp0); - // Normalized vector between two consecutive points - vector e10(p[f[fp1]] - p[f[fp0]]); + vector e10(p[f.nextLabel(fp0)] - p[f.thisLabel(fp0)]); scalar magE10 = mag(e10); e10 /= magE10 + VSMALL; diff --git a/src/parallel/decompose/faDecompose/faFieldDecomposer.C b/src/parallel/decompose/faDecompose/faFieldDecomposer.C index 6c6bc6a73b..59e7663f86 100644 --- a/src/parallel/decompose/faDecompose/faFieldDecomposer.C +++ b/src/parallel/decompose/faDecompose/faFieldDecomposer.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2016-2017 Wikki Ltd - Copyright (C) 2021 OpenCFD Ltd. + Copyright (C) 2021-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -126,7 +126,11 @@ processorAreaPatchFieldDecomposer mesh.edgeOwner(), mesh.edgeNeighbour(), addressingSlice, - mesh.weights().internalField() + ( + mesh.hasWeights() + ? mesh.weights().primitiveField() + : scalarField::null() + ) ) {} @@ -310,7 +314,9 @@ void Foam::faFieldDecomposer::reset(const faMesh& completeMesh) processorEdgePatchFieldDecomposerPtrs_.resize(nMappers); // Create weightings now - needed for proper parallel synchronization - (void)completeMesh.weights(); + //// (void)completeMesh.weights(); + // Disabled the above (2022-04-04) + // Use weights if they already exist, otherwise simply ignore // faPatches don't have their own start() - so these are invariant const labelList completePatchStarts