ENH: add proximityRegions filter to distanceSurface (#2108)

- combines region-based and proximity-based filtering

  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.

STYLE: 'proxityFaces' as newer synonym for 'proximity' filter

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

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<pointIndexHit>& 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<bool WantPointFilter = false>
 static bitSet simpleGeometricFilter
 (
     bitSet& ignoreCells,
     const List<pointIndexHit>& 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<false>(ignoreCells, nearest, fvmesh);
+            ignoreCellPoints = simpleGeometricFilter<false>
+            (
+                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<meshedSurface&>(*isoSurfacePtr_));
@@ -678,9 +721,13 @@ void Foam::distanceSurface::createGeometry()
         pointDistance_.clear();
     }
 
-    if (topologyFilterType::PROXIMITY == topoFilter_)
+    if (topologyFilterType::PROXIMITY_FACES == topoFilter_)
+    {
+        filterFaceProximity();
+    }
+    else if (topologyFilterType::PROXIMITY_REGIONS == topoFilter_)
     {
-        filterByProximity();
+        filterRegionProximity(ignoreCells);
     }
 
     if (debug)

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:
 

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<pointIndexHit> 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<scalar>());
+    Pstream::listCombineGather(areaRegion, plusEqOp<scalar>());
+
+    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<label>
+            (
+                meshCells_,
+                [&](const label celli){ return regionColour[celli]; }
+            )
+        );
+        writer.write("face-region", faceRegion);
+
+        // Region-wise filter state
+        labelField faceFilterState
+        (
+            ListOps::createWithValue<label>(surface_.size(), acceptFaces, 1, 0)
+        );
+        writer.write("filter-state", faceFilterState);
+    }
+
+
+    if (prune)
+    {
+        labelList pointMap, faceMap;
+        meshedSurface filtered
+        (
+            surface_.subsetMesh(acceptFaces, pointMap, faceMap)
+        );
+        surface_.transfer(filtered);
+
+        meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+    }
+}
+
 
+void Foam::distanceSurface::filterFaceProximity()
+{
+    const searchableSurface& geom = geometryPtr_();
+
+    // For face distances
+    const pointField& fc = surface_.faceCentres();
 
     // For each face
     scalarField faceDistance(fc.size(), GREAT);
@@ -356,38 +489,32 @@ void Foam::distanceSurface::filterByProximity()
 
     // Using the absolute proximity of the face centres is more robust.
 
+    bitSet acceptFaces(fc.size());
+    bool prune(false);
 
     // Consider the absolute proximity of the face centres
     forAll(faceDistance, facei)
     {
-        if (faceDistance[facei] <= absProximity_)
+        if (absProximity_ < faceDistance[facei])
         {
-            faceSelection.set(facei);
-        }
-        else
-        {
-            ++nTrimmed;
-
+            prune = true;
             if (debug & 2)
             {
                 Pout<< "trim reject: "
                     << faceDistance[facei] << nl;
             }
         }
+        else
+        {
+            acceptFaces.set(facei);
+        }
     }
 
 
     // Heavier debugging
     if (debug & 4)
     {
-        labelField faceFilterStatus(faceSelection.size(), Zero);
-
-        for (const label facei : faceSelection)
-        {
-            faceFilterStatus[facei] = 1;
-        }
-
-        const fileName outputName(surfaceName() + "-proximity-filter");
+        const fileName outputName(surfaceName() + "-face-proximity-filter");
 
         Info<< "Writing debug surface: " << outputName << nl;
 
@@ -400,16 +527,22 @@ void Foam::distanceSurface::filterByProximity()
 
         writer.write("absolute-distance", faceDistance);
         writer.write("normal-distance", faceNormalDistance);
-        writer.write("filter-state", faceFilterStatus);
+
+        // Region-wise filter state
+        labelField faceFilterState
+        (
+            ListOps::createWithValue<label>(surface_.size(), acceptFaces, 1, 0)
+        );
+        writer.write("filter-state", faceFilterState);
     }
 
 
-    if (returnReduce(nTrimmed, sumOp<label>()) != 0)
+    if (prune)
     {
         labelList pointMap, faceMap;
         meshedSurface filtered
         (
-            surface_.subsetMesh(faceSelection, pointMap, faceMap)
+            surface_.subsetMesh(acceptFaces, pointMap, faceMap)
         );
         surface_.transfer(filtered);