diff --git a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H index ce208cb259..8a365e9ae8 100644 --- a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H +++ b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2018-2020 OpenCFD Ltd. + Copyright (C) 2018-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -59,10 +59,10 @@ Usage bounds | Limit with bounding box | no | surfaceType | Type of surface | yes | surfaceName | Name of surface in \c triSurface/ | no | dict name - topology | Topology filter (none/largestRegion/nearestPoints/proximity) | no | none + topology | Topology filter name | no | none nearestPoints | Points for point-based segmentation | no | maxDistance | Max search distance for nearestPoints | no | GREAT - relProximity | Max limit of absolute vs normal distance | no | 1e3 + absProximity | Max proximity of face centres | no | 1e-5 \endtable SourceFiles diff --git a/src/sampling/surface/distanceSurface/distanceSurface.C b/src/sampling/surface/distanceSurface/distanceSurface.C index be57d3b122..3f634f04f9 100644 --- a/src/sampling/surface/distanceSurface/distanceSurface.C +++ b/src/sampling/surface/distanceSurface/distanceSurface.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2016-2020 OpenCFD Ltd. + Copyright (C) 2016-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -52,7 +52,9 @@ Foam::distanceSurface::topoFilterNames_ { topologyFilterType::NONE, "none" }, { topologyFilterType::LARGEST_REGION, "largestRegion" }, { topologyFilterType::NEAREST_POINTS, "nearestPoints" }, - { topologyFilterType::PROXIMITY, "proximity" }, + { topologyFilterType::PROXIMITY_REGIONS, "proximityRegions" }, + { topologyFilterType::PROXIMITY_FACES, "proximityFaces" }, + { topologyFilterType::PROXIMITY_FACES, "proximity" }, }); @@ -76,8 +78,9 @@ static inline void checkAllHits(const UList& nearest) if (notHit) { FatalErrorInFunction - << "Had " << notHit << " from " << nearest.size() - << " without a point hit" << endl + << "Had " << notHit << " faces/cells from " + << nearest.size() << " without a point hit." << nl + << "May be caused by a severely degenerate input surface" << nl << abort(FatalError); } } @@ -184,15 +187,18 @@ static inline void calcNormalDistance_filtered // the normal at the edge, which creates a zero-crossing // extending to infinity. // -// Ad hoc treatment: require that the surface hit -// point is within a somewhat generous bounding box -// for the cell +// Ad hoc treatment: require that the surface hit point is within a +// somewhat generous bounding box for the cell (+10%) +// +// Provisioning for filtering based on the cell points, +// but its usefulness remains to be seen (2020-12-09) template static bitSet simpleGeometricFilter ( bitSet& ignoreCells, const List& nearest, - const polyMesh& mesh + const polyMesh& mesh, + const scalar boundBoxInflate = 0.1 // 10% to catch corners ) { // A deny filter. Initially false (accept everything) @@ -215,9 +221,7 @@ static bitSet simpleGeometricFilter cellBb.clear(); cellBb.add(mesh.points(), cPoints); - - // Expand slightly to catch corners - cellBb.inflate(0.1); + cellBb.inflate(boundBoxInflate); if (!cellBb.contains(pt)) { @@ -455,10 +459,30 @@ void Foam::distanceSurface::createGeometry() checkAllHits(nearest); // Geometric pre-filtering when distance == 0 + + // NOTE (2021-05-31) + // Can skip the prefilter if we use proximity-regions filter anyhow + // but it makes the iso algorithm more expensive and doesn't help + // unless we start relying on area-based weighting for rejecting regions. + if (filterCells) { - ignoreCellPoints = - simpleGeometricFilter(ignoreCells, nearest, fvmesh); + ignoreCellPoints = simpleGeometricFilter + ( + ignoreCells, + nearest, + fvmesh, + + // Inflate bound box. + // - To catch corners: approx. 10% + // - Extra generous for PROXIMITY_REGIONS + // (extra weighting for 'bad' faces) + ( + topologyFilterType::PROXIMITY_REGIONS == topoFilter_ + ? 1 + : 0.1 + ) + ); } if (withSignDistance_) @@ -617,6 +641,8 @@ void Foam::distanceSurface::createGeometry() refineBlockedCells(ignoreCells, isoCutter); filterKeepNearestRegions(ignoreCells); } + + // Note: apply similar filtering for PROXIMITY_REGIONS later instead } // Don't need point filter beyond this point @@ -627,6 +653,7 @@ void Foam::distanceSurface::createGeometry() { Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl; cellDistance.write(); + pointScalarField pDist ( IOobject @@ -660,6 +687,21 @@ void Foam::distanceSurface::createGeometry() ); + // Restrict ignored cells to those actually cut + if (filterCells && topologyFilterType::PROXIMITY_REGIONS == topoFilter_) + { + isoSurfaceBase isoCutter + ( + mesh_, + cellDistance, + pointDistance_, + distance_ + ); + + refineBlockedCells(ignoreCells, isoCutter); + } + + // ALGO_POINT still needs cell, point fields (for interpolate) // The others can do straight transfer @@ -667,7 +709,8 @@ void Foam::distanceSurface::createGeometry() if ( isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT - || topologyFilterType::PROXIMITY == topoFilter_ + || topologyFilterType::PROXIMITY_FACES == topoFilter_ + || topologyFilterType::PROXIMITY_REGIONS == topoFilter_ ) { surface_.transfer(static_cast(*isoSurfacePtr_)); @@ -678,9 +721,13 @@ void Foam::distanceSurface::createGeometry() pointDistance_.clear(); } - if (topologyFilterType::PROXIMITY == topoFilter_) + if (topologyFilterType::PROXIMITY_FACES == topoFilter_) { - filterByProximity(); + filterFaceProximity(); + } + else if (topologyFilterType::PROXIMITY_REGIONS == topoFilter_) + { + filterRegionProximity(ignoreCells); } if (debug) diff --git a/src/sampling/surface/distanceSurface/distanceSurface.H b/src/sampling/surface/distanceSurface/distanceSurface.H index 2273839329..621ad8490b 100644 --- a/src/sampling/surface/distanceSurface/distanceSurface.H +++ b/src/sampling/surface/distanceSurface/distanceSurface.H @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2016 OpenFOAM Foundation - Copyright (C) 2016-2020 OpenCFD Ltd. + Copyright (C) 2016-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -45,7 +45,7 @@ Usage type distanceSurface; surfaceType triSurfaceMesh; surfaceName something.obj; - topology proximity; + topology proximityFaces; } surface2 @@ -77,27 +77,37 @@ Usage bounds | Limit with bounding box | no | surfaceType | Type of surface | yes | surfaceName | Name of surface in \c triSurface/ | no | dict name - topology | Topology filter (none/largestRegion/nearestPoints/proximity) | no | none + topology | Topology filter name | no | none nearestPoints | Points for point-based segmentation | no | maxDistance | Max search distance for nearestPoints | no | GREAT absProximity | Max proximity of face centres | no | 1e-5 \endtable - Filtering (for zero-distance only) + Topology/Filtering (for zero-distance only). + These represent different ways to tackle the "ragged edge" problem. + + - \c none : + No filtering + + - \c proximityFaces or \c proximity (post-filter): + Checks the resulting faces against the original search surface + and rejects faces with a distance greater than \c absProximity. + + - \c proxityRegions (post-filter): + Checks the distance of the resulting faces against the original + search surface. Filters based on the area-weighted distance + of each topologically connected region. + If the area-weighted distance of a region is greater than + \c absProximity, the entire region is rejected. - \c largestRegion (pre-filter): The cut cells are checked for topological connectivity and the region with the most number of cut cells is retained. - This handles the "ragged" edge problem. - \c nearestPoints (pre-filter): The cut cells split into regions, the regions closest to the user-defined points are retained. Uses \c maxDistance for additional control. - - - \c proximity (post-filter): - Checks the resulting faces against the original search surface - and rejects faces with a distance greater than \c absProximity. . Note @@ -116,6 +126,7 @@ Changed default algorithm from cell to topo (2020-12). SourceFiles distanceSurface.C + distanceSurfaceFilter.C \*---------------------------------------------------------------------------*/ @@ -146,7 +157,9 @@ class distanceSurface NONE, //!< No additional filtering LARGEST_REGION, //!< Retain largest region NEAREST_POINTS, //!< Retain regions nearest to the points - PROXIMITY //!< Post-filter for surface proximity + PROXIMITY_REGIONS, //!< Retain regions with good surface proximity + PROXIMITY_FACES, //!< Retain faces with good surface proximity + PROXIMITY = PROXIMITY_FACES }; //- Names for the topology filter @@ -262,15 +275,18 @@ protected: //- Prepare blockedFaces for region split bitSet filterPrepareRegionSplit(const bitSet& ignoreCells) const; - //- Adjust extracted iso-surface to remove far faces - void filterByProximity(); - //- Keep region with the most cuts (after region split) void filterKeepLargestRegion(bitSet& ignoreCells) const; //- Keep region(s) closest to the nearest points void filterKeepNearestRegions(bitSet& ignoreCells) const; + //- Remove region(s) that have far faces + void filterRegionProximity(bitSet& ignoreCells) const; + + //- Adjust extracted iso-surface to remove far faces + void filterFaceProximity(); + public: diff --git a/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C b/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C index 9599058302..0ca4e4961b 100644 --- a/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C +++ b/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2020 OpenCFD Ltd. + Copyright (C) 2020-2021 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,6 +28,7 @@ License #include "distanceSurface.H" #include "regionSplit.H" #include "syncTools.H" +#include "ListOps.H" #include "vtkSurfaceWriter.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -306,16 +307,148 @@ void Foam::distanceSurface::filterKeepNearestRegions } -void Foam::distanceSurface::filterByProximity() +void Foam::distanceSurface::filterRegionProximity +( + bitSet& ignoreCells +) const { const searchableSurface& geom = geometryPtr_(); - // Filtering for faces + // For face distances const pointField& fc = surface_.faceCentres(); - bitSet faceSelection(fc.size()); - label nTrimmed = 0; + // For region split + bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells)); + // Split region + regionSplit rs(mesh_, blockedFaces); + blockedFaces.clearStorage(); + + const labelList& regionColour = rs; + + // For each face + scalarField faceDistance(fc.size(), GREAT); + + { + List nearest; + geom.findNearest + ( + fc, + // Use initialized field (GREAT) to limit search too + faceDistance, + nearest + ); + calcAbsoluteDistance(faceDistance, fc, nearest); + } + + // Identical number of regions on all processors + scalarField areaRegion(rs.nRegions(), Zero); + scalarField distRegion(rs.nRegions(), Zero); + + forAll(meshCells_, facei) + { + const label celli = meshCells_[facei]; + const label regioni = regionColour[celli]; + + const scalar faceArea = surface_[facei].mag(surface_.points()); + distRegion[regioni] += (faceDistance[facei] * faceArea); + areaRegion[regioni] += (faceArea); + } + + Pstream::listCombineGather(distRegion, plusEqOp()); + Pstream::listCombineGather(areaRegion, plusEqOp()); + + if (Pstream::master()) + { + forAll(distRegion, regioni) + { + distRegion[regioni] /= (areaRegion[regioni] + VSMALL); + } + } + + Pstream::listCombineScatter(distRegion); + + + // Define the per-face acceptance based on the region average distance + + bitSet acceptFaces(fc.size()); + bool prune(false); + + forAll(meshCells_, facei) + { + const label celli = meshCells_[facei]; + const label regioni = regionColour[celli]; + + // NB: do not filter by individual faces as well since this + // has been reported to cause minor holes for surfaces with + // high curvature! (2021-06-10) + + if (absProximity_ < distRegion[regioni]) + { + prune = true; + } + else + { + acceptFaces.set(facei); + } + } + + // Heavier debugging + if (debug & 4) + { + const fileName outputName(surfaceName() + "-region-proximity-filter"); + + Info<< "Writing debug surface: " << outputName << nl; + + surfaceWriters::vtkWriter writer + ( + surface_.points(), + surface_, // faces + outputName + ); + + writer.write("absolute-distance", faceDistance); + + // Region segmentation + labelField faceRegion + ( + ListOps::create