From 18b134319d93c635a7765811ec3bfd263b38ebd8 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Mon, 23 Jul 2018 15:18:20 +0200 Subject: [PATCH] ENH: handling of open edges distanceSurface at 0 (issue #948) Some special adjustments are undertaken for distance = 0. - With the isoSurfaceCell algorithm is used, additional checks for open surfaces edges are used to limit the extend of resulting distance surface. The resulting surface elements will not, however, contain partial cell coverage. - Always treated as signed (ignoring the input value), since it is nearly impossible to generate any surface otherwise. --- .../surface/distanceSurface/distanceSurface.C | 59 +++++++++++++-- .../surface/distanceSurface/distanceSurface.H | 10 ++- src/sampling/surface/isoSurface/isoSurface.C | 73 +++++++++---------- src/sampling/surface/isoSurface/isoSurface.H | 9 ++- .../surface/isoSurface/isoSurfaceCell.C | 36 +++++---- .../surface/isoSurface/isoSurfaceCell.H | 24 ++++-- 6 files changed, 139 insertions(+), 72 deletions(-) diff --git a/src/sampling/surface/distanceSurface/distanceSurface.C b/src/sampling/surface/distanceSurface/distanceSurface.C index e22cf4fd8b..b49b191773 100644 --- a/src/sampling/surface/distanceSurface/distanceSurface.C +++ b/src/sampling/surface/distanceSurface/distanceSurface.C @@ -66,7 +66,11 @@ Foam::distanceSurface::distanceSurface ) ), distance_(dict.get("distance")), - signed_(dict.get("signed")), + signed_ + ( + // Always signed when distance = 0. + dict.get("signed") || equal(distance_, Zero) + ), cell_(dict.lookupOrDefault("cell", true)), regularise_(dict.lookupOrDefault("regularise", true)), bounds_(dict.lookupOrDefault("bounds", boundBox::invertedBox)), @@ -108,7 +112,11 @@ Foam::distanceSurface::distanceSurface ) ), distance_(distance), - signed_(signedDistance), + signed_ + ( + // Always signed when distance = 0. + signedDistance || equal(distance_, Zero) + ), cell_(cell), regularise_(regularise), bounds_(bounds), @@ -155,6 +163,17 @@ void Foam::distanceSurface::createGeometry() ); volScalarField& cellDistance = *cellDistancePtr_; + // For distance = 0 (and isoSurfaceCell) we apply additional filtering + // to limit the extent of open edges. + + const bool filterCells = (cell_ && signed_ && equal(distance_, Zero)); + + bitSet ignoreCells; + if (filterCells) + { + ignoreCells.resize(fvm.nCells()); + } + // Internal field { const pointField& cc = fvm.C(); @@ -173,11 +192,26 @@ void Foam::distanceSurface::createGeometry() vectorField norms; surfPtr_().getNormal(nearest, norms); + boundBox cellBb; + forAll(norms, i) { const point diff(cc[i] - nearest[i].hitPoint()); fld[i] = sign(diff & norms[i]) * Foam::mag(diff); + + // Since cellPoints() is used in isoSurfaceCell too, + // no additional overhead caused here + if (filterCells) + { + cellBb.clear(); + cellBb.add(fvm.points(), fvm.cellPoints(i)); + + if (!cellBb.contains(nearest[i].hitPoint())) + { + ignoreCells.set(i); + } + } } } else @@ -185,15 +219,17 @@ void Foam::distanceSurface::createGeometry() forAll(nearest, i) { fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint()); + + // No filtering for unsigned or distance != 0. } } } - volScalarField::Boundary& cellDistanceBf = - cellDistance.boundaryFieldRef(); - // Patch fields { + volScalarField::Boundary& cellDistanceBf = + cellDistance.boundaryFieldRef(); + forAll(fvm.C().boundaryField(), patchi) { const pointField& cc = fvm.C().boundaryField()[patchi]; @@ -293,6 +329,12 @@ void Foam::distanceSurface::createGeometry() pDist.write(); } + // Don't need ignoreCells if there is nothing to ignore. + if (!ignoreCells.any()) + { + ignoreCells.clear(); + } + // Direct from cell field and point field. if (cell_) @@ -306,7 +348,9 @@ void Foam::distanceSurface::createGeometry() pointDistance_, distance_, regularise_, - bounds_ + bounds_, + 1e-6, // mergeTol + ignoreCells ) ); } @@ -320,7 +364,8 @@ void Foam::distanceSurface::createGeometry() pointDistance_, distance_, regularise_, - bounds_ + bounds_, + 1e-6 ) ); } diff --git a/src/sampling/surface/distanceSurface/distanceSurface.H b/src/sampling/surface/distanceSurface/distanceSurface.H index c20166e6dd..e0f59de1f6 100644 --- a/src/sampling/surface/distanceSurface/distanceSurface.H +++ b/src/sampling/surface/distanceSurface/distanceSurface.H @@ -38,9 +38,17 @@ Usage regularise | point snapping | yes | bounds | limit with bounding box | no | surfaceType | type of surface | yes | - surfaceName | name of surface in triSurface/ | no | dict name + surfaceName | name of surface in \c triSurface/ | no | dict name \endtable +Note + Some special adjustments are undertaken for distance = 0. + - Always treated as signed (ignoring the input value), since it is + nearly impossible to generate any surface otherwise. + - With the isoSurfaceCell algorithm is used, additional checks for open + surfaces edges are used to limit the extend of resulting distance + surface. The resulting surface elements will not, however, contain + partial cell coverage. SourceFiles distanceSurface.C diff --git a/src/sampling/surface/isoSurface/isoSurface.C b/src/sampling/surface/isoSurface/isoSurface.C index d3faf431f2..c5ae60b669 100644 --- a/src/sampling/surface/isoSurface/isoSurface.C +++ b/src/sampling/surface/isoSurface/isoSurface.C @@ -138,17 +138,17 @@ void Foam::isoSurface::syncUnseparatedPoints if (Pstream::parRun()) { // Send - forAll(patches, patchi) + for (const polyPatch& p : patches) { if ( - isA(patches[patchi]) - && patches[patchi].nPoints() > 0 - && collocatedPatch(patches[patchi]) + isA(p) + && p.nPoints() > 0 + && collocatedPatch(p) ) { const processorPolyPatch& pp = - refCast(patches[patchi]); + refCast(p); const labelList& meshPts = pp.meshPoints(); const labelList& nbrPts = pp.neighbPoints(); @@ -171,18 +171,17 @@ void Foam::isoSurface::syncUnseparatedPoints } // Receive and combine. - - forAll(patches, patchi) + for (const polyPatch& p : patches) { if ( - isA(patches[patchi]) - && patches[patchi].nPoints() > 0 - && collocatedPatch(patches[patchi]) + isA(p) + && p.nPoints() > 0 + && collocatedPatch(p) ) { const processorPolyPatch& pp = - refCast(patches[patchi]); + refCast(p); pointField nbrPatchInfo(pp.nPoints()); { @@ -200,7 +199,7 @@ void Foam::isoSurface::syncUnseparatedPoints forAll(meshPts, pointi) { - label meshPointi = meshPts[pointi]; + const label meshPointi = meshPts[pointi]; minEqOp() ( pointValues[meshPointi], @@ -212,12 +211,12 @@ void Foam::isoSurface::syncUnseparatedPoints } // Do the cyclics. - forAll(patches, patchi) + for (const polyPatch& p : patches) { - if (isA(patches[patchi])) + if (isA(p)) { const cyclicPolyPatch& cycPatch = - refCast(patches[patchi]); + refCast(p); if (cycPatch.owner() && collocatedPatch(cycPatch)) { @@ -241,8 +240,8 @@ void Foam::isoSurface::syncUnseparatedPoints forAll(coupledPoints, i) { const edge& e = coupledPoints[i]; - label p0 = meshPts[e[0]]; - label p1 = nbrMeshPoints[e[1]]; + const label p0 = meshPts[e[0]]; + const label p1 = nbrMeshPoints[e[1]]; minEqOp()(pointValues[p0], half1Values[i]); minEqOp()(pointValues[p1], half0Values[i]); @@ -410,10 +409,8 @@ void Foam::isoSurface::calcCutTypes } } - forAll(patches, patchi) + for (const polyPatch& pp : patches) { - const polyPatch& pp = patches[patchi]; - label facei = pp.start(); forAll(pp, i) @@ -1339,8 +1336,8 @@ Foam::triSurface Foam::isoSurface::subsetMesh Foam::isoSurface::isoSurface ( - const volScalarField& cVals, - const scalarField& pVals, + const volScalarField& cellValues, + const scalarField& pointValues, const scalar iso, const bool regularise, const boundBox& bounds, @@ -1348,8 +1345,8 @@ Foam::isoSurface::isoSurface ) : MeshStorage(), - mesh_(cVals.mesh()), - pVals_(pVals), + mesh_(cellValues.mesh()), + pVals_(pointValues), iso_(iso), regularise_(regularise), bounds_(bounds), @@ -1358,10 +1355,10 @@ Foam::isoSurface::isoSurface if (debug) { Pout<< "isoSurface:" << nl - << " isoField : " << cVals.name() << nl + << " isoField : " << cellValues.name() << nl << " cell min/max : " - << min(cVals.primitiveField()) << " / " - << max(cVals.primitiveField()) << nl + << min(cellValues.primitiveField()) << " / " + << max(cellValues.primitiveField()) << nl << " point min/max : " << min(pVals_) << " / " << max(pVals_) << nl @@ -1379,7 +1376,7 @@ Foam::isoSurface::isoSurface // Rewrite input volScalarField to have interpolated values // on separated patches. - cValsPtr_.reset(adaptPatchFields(cVals).ptr()); + cValsPtr_.reset(adaptPatchFields(cellValues).ptr()); // Construct cell centres field consistent with cVals @@ -1411,7 +1408,7 @@ Foam::isoSurface::isoSurface // Adapt separated coupled (proc and cyclic) patches if (pp.coupled()) { - fvPatchVectorField& pfld = const_cast + auto& pfld = const_cast ( meshC.boundaryField()[patchi] ); @@ -1431,9 +1428,10 @@ Foam::isoSurface::isoSurface } else if (isA(pp)) { - typedef slicedVolVectorField::Boundary bType; - - bType& bfld = const_cast(meshC.boundaryField()); + auto& bfld = const_cast + ( + meshC.boundaryField() + ); // Clear old value. Cannot resize it since is a slice. bfld.set(patchi, nullptr); @@ -1551,18 +1549,15 @@ Foam::isoSurface::isoSurface // Determine if point is on boundary. bitSet isBoundaryPoint(mesh_.nPoints()); - forAll(patches, patchi) + for (const polyPatch& pp : patches) { // Mark all boundary points that are not physically coupled // (so anything but collocated coupled patches) - if (patches[patchi].coupled()) + if (pp.coupled()) { const coupledPolyPatch& cpp = - refCast - ( - patches[patchi] - ); + refCast(pp); bitSet isCollocated(collocatedFaces(cpp)); @@ -1578,8 +1573,6 @@ Foam::isoSurface::isoSurface } else { - const polyPatch& pp = patches[patchi]; - forAll(pp, i) { const face& f = mesh_.faces()[pp.start()+i]; diff --git a/src/sampling/surface/isoSurface/isoSurface.H b/src/sampling/surface/isoSurface/isoSurface.H index 1fd72e6550..3df599257e 100644 --- a/src/sampling/surface/isoSurface/isoSurface.H +++ b/src/sampling/surface/isoSurface/isoSurface.H @@ -406,14 +406,17 @@ public: //- Construct from cell values and point values. // Uses boundaryField for boundary values. // Holds reference to cellIsoVals and pointIsoVals. + // + // \param bounds optional bounding box for trimming + // \param mergeTol fraction of mesh bounding box for merging points isoSurface ( - const volScalarField& cellIsoVals, - const scalarField& pointIsoVals, + const volScalarField& cellValues, + const scalarField& pointValues, const scalar iso, const bool regularise, const boundBox& bounds = boundBox::invertedBox, - const scalar mergeTol = 1e-6 // fraction of bounding box + const scalar mergeTol = 1e-6 ); diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.C b/src/sampling/surface/isoSurface/isoSurfaceCell.C index 47b539f250..f63c5e7c1c 100644 --- a/src/sampling/surface/isoSurface/isoSurfaceCell.C +++ b/src/sampling/surface/isoSurface/isoSurfaceCell.C @@ -85,6 +85,11 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType const label celli ) const { + if (ignoreCells_.test(celli)) + { + return NOTCUT; + } + const cell& cFaces = mesh_.cells()[celli]; if (isTet.test(celli)) @@ -174,8 +179,9 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType { return SPHERE; } - else + else if (nCuts > 1) { + // Need at least two edge cuts, otherwise this is a spurious cut return CUT; } } @@ -1268,20 +1274,22 @@ Foam::triSurface Foam::isoSurfaceCell::subsetMesh Foam::isoSurfaceCell::isoSurfaceCell ( const polyMesh& mesh, - const scalarField& cVals, - const scalarField& pVals, + const scalarField& cellValues, + const scalarField& pointValues, const scalar iso, const bool regularise, const boundBox& bounds, - const scalar mergeTol + const scalar mergeTol, + const bitSet& ignoreCells ) : MeshStorage(), mesh_(mesh), - cVals_(cVals), - pVals_(pVals), + cVals_(cellValues), + pVals_(pointValues), iso_(iso), bounds_(bounds), + ignoreCells_(ignoreCells), mergeDistance_(mergeTol*mesh.bounds().mag()) { if (debug) @@ -1298,6 +1306,8 @@ Foam::isoSurfaceCell::isoSurfaceCell << " mergeTol : " << mergeTol << nl << " mesh span : " << mesh.bounds().mag() << nl << " mergeDistance : " << mergeDistance_ << nl + << " ignoreCells : " << ignoreCells_.count() + << " / " << cVals_.size() << nl << endl; } @@ -1317,7 +1327,7 @@ Foam::isoSurfaceCell::isoSurfaceCell // Determine if any cut through cell - calcCutTypes(isTet, cVals, pVals); + calcCutTypes(isTet, cellValues, pointValues); if (debug && isA(mesh)) { @@ -1361,8 +1371,8 @@ Foam::isoSurfaceCell::isoSurfaceCell calcSnappedCc ( isTet, - cVals, - pVals, + cellValues, + pointValues, snappedPoints, snappedCc ); @@ -1389,8 +1399,8 @@ Foam::isoSurfaceCell::isoSurfaceCell calcSnappedPoint ( isTet, - cVals, - pVals, + cellValues, + pointValues, snappedPoints, snappedPoint ); @@ -1417,8 +1427,8 @@ Foam::isoSurfaceCell::isoSurfaceCell generateTriPoints ( - cVals, - pVals, + cellValues, + pointValues, mesh_.cellCentres(), mesh_.points(), diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.H b/src/sampling/surface/isoSurface/isoSurfaceCell.H index 57b4d1e743..3738253c3a 100644 --- a/src/sampling/surface/isoSurface/isoSurfaceCell.H +++ b/src/sampling/surface/isoSurface/isoSurfaceCell.H @@ -78,16 +78,16 @@ class isoSurfaceCell enum segmentCutType { - NEARFIRST, // intersection close to e.first() - NEARSECOND, // ,, e.second() - NOTNEAR // no intersection + NEARFIRST, //!< Intersection close to e.first() + NEARSECOND, //!< Intersection close to e.second() + NOTNEAR //!< No intersection }; enum cellCutType { - NOTCUT, // not cut - SPHERE, // all edges to cell centre cut - CUT // normal cut + NOTCUT, //!< Cell not cut + SPHERE, //!< All edges to cell centre cut + CUT //!< Normal cut }; @@ -104,6 +104,9 @@ class isoSurfaceCell //- Optional bounds const boundBox bounds_; + //- Optional cells to ignore. + const bitSet& ignoreCells_; + //- When to merge points const scalar mergeDistance_; @@ -324,7 +327,11 @@ public: // Constructors - //- Construct from dictionary + //- Construct from cell and point values + // + // \param bounds optional bounding box for trimming + // \param mergeTol fraction of mesh bounding box for merging points + // \param ignoreCells cells to ignore in the cellValues isoSurfaceCell ( const polyMesh& mesh, @@ -333,7 +340,8 @@ public: const scalar iso, const bool regularise, const boundBox& bounds = boundBox::invertedBox, - const scalar mergeTol = 1e-6 // fraction of bounding box + const scalar mergeTol = 1e-6, + const bitSet& ignoreCells = bitSet() );