zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-distanceSurface-filtering' into 'develop'

Browse Source 
ENH: add proximityRegions filter to distanceSurface (#2108)

See merge request Development/openfoam!460
This commit is contained in:
Andrew Heather
2021-06-10 19:13:40 +00:00
parent b6b4ab071d ca88abba7e
commit 596d8fab39
4 changed files with 251 additions and 55 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
View File

@ -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

79
src/sampling/surface/distanceSurface/distanceSurface.C
View File

@ -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_)
{
filterByProximity();
filterFaceProximity();
}
else if (topologyFilterType::PROXIMITY_REGIONS == topoFilter_)
{
filterRegionProximity(ignoreCells);
}
if (debug)

42
src/sampling/surface/distanceSurface/distanceSurface.H
View File

@ -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:

179
src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
View File

@ -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);