From 2f6082712e7450a0433de6a4900b4ff19a1f3839 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Mon, 30 Nov 2020 13:05:56 +0100 Subject: [PATCH] ENH: modernize some code constructs in isoSurface - add debug field to isoSurfaceTopo - don't need dynamic field for new points - reduce code in sampledIsoSurfaceCell --- .../isoSurface/sampledIsoSurfaceCell.C | 9 +- .../surface/isoSurface/isoSurfaceBase.H | 26 ++ .../isoSurface/isoSurfaceCellTemplates.C | 2 +- .../surface/isoSurface/isoSurfacePoint.C | 165 ++++----- .../surface/isoSurface/isoSurfacePoint.H | 5 +- .../isoSurface/isoSurfacePointTemplates.C | 2 +- .../surface/isoSurface/isoSurfaceTopo.C | 313 ++++++++---------- .../surface/isoSurface/isoSurfaceTopo.H | 22 +- .../isoSurface/isoSurfaceTopoTemplates.C | 41 ++- 9 files changed, 283 insertions(+), 302 deletions(-) diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C index fbade62037..4a5c8ff97b 100644 --- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C +++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C @@ -139,6 +139,9 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const } } + + meshedSurface& mySurface = const_cast(*this); + isoSurfaceCell surf ( fvm, @@ -148,11 +151,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const isoParams_ ); - // Replace current geomety - const_cast - ( - *this - ).transfer(static_cast(surf)); + mySurface.transfer(static_cast(surf)); meshCells_.transfer(surf.meshCells()); if (debug) diff --git a/src/sampling/surface/isoSurface/isoSurfaceBase.H b/src/sampling/surface/isoSurface/isoSurfaceBase.H index e225912de1..3229f002f5 100644 --- a/src/sampling/surface/isoSurface/isoSurfaceBase.H +++ b/src/sampling/surface/isoSurface/isoSurfaceBase.H @@ -54,6 +54,9 @@ SourceFiles namespace Foam { +// Forward Declarations +class tetCell; + /*---------------------------------------------------------------------------*\ Class isoSurfaceBase Declaration \*---------------------------------------------------------------------------*/ @@ -77,6 +80,29 @@ protected: labelList meshCells_; + // Protected Member Functions + + //- Check for tet values above/below given (iso) value + // Result encoded as a single integer + inline static constexpr int getTetCutIndex + ( + const scalar a, + const scalar b, + const scalar c, + const scalar d, + const scalar isoval + ) noexcept + { + return + ( + (a < isoval ? 1 : 0) + | (b < isoval ? 2 : 0) + | (c < isoval ? 4 : 0) + | (d < isoval ? 8 : 0) + ); + } + + public: // Constructors diff --git a/src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C index 5402b396c8..71f033cd76 100644 --- a/src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C +++ b/src/sampling/surface/isoSurface/isoSurfaceCellTemplates.C @@ -68,7 +68,7 @@ Type Foam::isoSurfaceCell::generatePoint } else { - scalar s = 0.4999; + constexpr scalar s = 0.4999; return s*p1 + (1.0-s)*p0; } diff --git a/src/sampling/surface/isoSurface/isoSurfacePoint.C b/src/sampling/surface/isoSurface/isoSurfacePoint.C index 73bd9d7de9..9cc8eada00 100644 --- a/src/sampling/surface/isoSurface/isoSurfacePoint.C +++ b/src/sampling/surface/isoSurface/isoSurfacePoint.C @@ -46,22 +46,53 @@ namespace Foam namespace Foam { - // Helper class for slicing triangles - class storeOp + +// Helper class for slicing triangles +struct storeOp +{ + DynamicList& list; + + storeOp(DynamicList& tris) + : + list(tris) + {} + + void operator()(const triPoints& tri) { - public: - DynamicList& tris_; + list.append(tri); + } - inline storeOp(DynamicList& tris) - : - tris_(tris) - {} + void operator()(triPoints&& tri) + { + list.append(std::move(tri)); + } +}; - inline void operator()(const triPoints& tri) - { - tris_.append(tri); - } - }; + +// Is patch co-located (i.e. non-separated) coupled patch? +static inline bool collocatedPatch(const polyPatch& pp) +{ + const auto* cpp = isA(pp); + return bool(cpp) && cpp->parallel() && !cpp->separated(); +} + + +// Collocated patch, with extra checks +#if 0 +static bool isCollocatedPatch(const coupledPolyPatch& pp) +{ + if (isA(pp) || isA(pp)) + { + return collocatedPatch(pp); + } + + FatalErrorInFunction + << "Unhandled coupledPolyPatch type " << pp.type() << nl + << abort(FatalError); + + return false; +} +#endif } // End namespace Foam @@ -83,39 +114,20 @@ bool Foam::isoSurfacePoint::noTransform(const tensor& tt) const } -bool Foam::isoSurfacePoint::collocatedPatch(const polyPatch& pp) -{ - const coupledPolyPatch& cpp = refCast(pp); - return cpp.parallel() && !cpp.separated(); -} - - Foam::bitSet Foam::isoSurfacePoint::collocatedFaces ( const coupledPolyPatch& pp -) const +) { // Initialise to false bitSet collocated(pp.size()); - if (isA(pp)) + if (isA(pp) || isA(pp)) { if (collocatedPatch(pp)) { - forAll(pp, i) - { - collocated.set(i); - } - } - } - else if (isA(pp)) - { - if (collocatedPatch(pp)) - { - forAll(pp, i) - { - collocated.set(i); - } + // All on + collocated = true; } } else @@ -377,7 +389,7 @@ void Foam::isoSurfacePoint::calcCutTypes for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei) { - const scalar& ownValue = cVals[own[facei]]; + const scalar ownValue = cVals[own[facei]]; const bool ownLower = (ownValue < iso_); scalar nbrValue; @@ -403,7 +415,7 @@ void Foam::isoSurfacePoint::calcCutTypes { // Is mesh edge cut? // - by looping over all the edges of the face. - const face f = mesh_.faces()[facei]; + const face& f = mesh_.faces()[facei]; if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower)) { @@ -414,11 +426,9 @@ void Foam::isoSurfacePoint::calcCutTypes for (const polyPatch& pp : patches) { - label facei = pp.start(); - - forAll(pp, i) + for (const label facei : pp.range()) { - const scalar& ownValue = cVals[own[facei]]; + const scalar ownValue = cVals[own[facei]]; const bool ownLower = (ownValue < iso_); scalar nbrValue; @@ -444,15 +454,13 @@ void Foam::isoSurfacePoint::calcCutTypes { // Is mesh edge cut? // - by looping over all the edges of the face. - const face f = mesh_.faces()[facei]; + const face& f = mesh_.faces()[facei]; if (isEdgeOfFaceCut(pVals, f, ownLower, neiLower)) { faceCutType_[facei] = CUT; } } - - ++facei; } } @@ -545,9 +553,7 @@ void Foam::isoSurfacePoint::calcSnappedCc { if (cellCutType_[celli] == CUT) { - scalar cVal = cVals[celli]; - - const cell& cFaces = mesh_.cells()[celli]; + const scalar cVal = cVals[celli]; localTriPoints.clear(); label nOther = 0; @@ -556,9 +562,9 @@ void Foam::isoSurfacePoint::calcSnappedCc // Create points for all intersections close to cell centre // (i.e. from pyramid edges) - forAll(cFaces, cFacei) + for (const label facei : mesh_.cells()[celli]) { - label facei = cFaces[cFacei]; + const face& f = mesh_.faces()[facei]; scalar nbrValue; point nbrPoint; @@ -581,8 +587,6 @@ void Foam::isoSurfacePoint::calcSnappedCc s[2] = isoFraction(cVal, nbrValue); pt[2] = (1.0-s[2])*cc[celli] + s[2]*nbrPoint; - const face& f = mesh_.faces()[cFaces[cFacei]]; - forAll(f, fp) { // From cc to fp @@ -638,10 +642,9 @@ void Foam::isoSurfacePoint::calcSnappedCc else if (localTriPoints.size() == 3) { // Single triangle. No need for any analysis. Average points. - pointField points; - points.transfer(localTriPoints); snappedCc[celli] = snappedPoints.size(); - snappedPoints.append(sum(points)/points.size()); + snappedPoints.append(sum(localTriPoints)/3); + localTriPoints.clear(); //Pout<< " point:" << pointi // << " replacing coord:" << mesh_.points()[pointi] @@ -741,7 +744,7 @@ void Foam::isoSurfacePoint::calcSnappedPoint // Create points for all intersections close to point // (i.e. from pyramid edges) const face& f = mesh_.faces()[facei]; - label own = mesh_.faceOwner()[facei]; + const label own = mesh_.faceOwner()[facei]; // Get neighbour value scalar nbrValue; @@ -1087,19 +1090,17 @@ void Foam::isoSurfacePoint::trimToPlanes if (useA) { - insideOpB.tris_.clear(); - forAll(insideOpA.tris_, i) + insideTrisB.clear(); + for (const triPoints& tri : insideTrisA) { - const triPoints& tri = insideOpA.tris_[i]; triPointRef(tri).sliceWithPlane(pln, insideOpB, dop); } } else { - insideOpA.tris_.clear(); - forAll(insideOpB.tris_, i) + insideTrisA.clear(); + for (const triPoints& tri : insideTrisB) { - const triPoints& tri = insideOpB.tris_[i]; triPointRef(tri).sliceWithPlane(pln, insideOpA, dop); } } @@ -1107,26 +1108,16 @@ void Foam::isoSurfacePoint::trimToPlanes } + DynamicList& newTris = (useA ? insideTrisA : insideTrisB); + + newTriPoints.reserve(newTriPoints.size() + 3*newTris.size()); + // Transfer - if (useA) + for (const triPoints& tri : newTris) { - forAll(insideOpA.tris_, i) - { - const triPoints& tri = insideOpA.tris_[i]; - newTriPoints.append(tri[0]); - newTriPoints.append(tri[1]); - newTriPoints.append(tri[2]); - } - } - else - { - forAll(insideOpB.tris_, i) - { - const triPoints& tri = insideOpB.tris_[i]; - newTriPoints.append(tri[0]); - newTriPoints.append(tri[1]); - newTriPoints.append(tri[2]); - } + newTriPoints.append(tri[0]); + newTriPoints.append(tri[1]); + newTriPoints.append(tri[2]); } } @@ -1449,25 +1440,15 @@ Foam::isoSurfacePoint::isoSurfacePoint } - // Pre-calculate patch-per-face to avoid whichPatch call. labelList boundaryRegion(mesh_.nBoundaryFaces()); - forAll(patches, patchi) + for (const polyPatch& pp : patches) { - const polyPatch& pp = patches[patchi]; - - label facei = pp.start(); - - forAll(pp, i) - { - boundaryRegion[facei-mesh_.nInternalFaces()] = patchi; - ++facei; - } + SubList