From 4c3a95c808bc0970ebb5ae38232687fe8014ad27 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Wed, 9 Dec 2020 11:36:49 +0100
Subject: [PATCH] ENH: additional filtering for distance surface (#1950)
- adds topology-based segmentation of the surfaces generated with
distance surfaces. This can occur when the surface terminates
close to a thin wall gap in the mesh; resulting in a cuts that
extend into the next region.
The cutting algorithm does not normally distinguish between these
types of "ragged" cuts, and legitimate ones (eg, cutting multiple
pipes). The additional segmentation controls provide for two common
scenarios:
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.
nearestPoints (pre-filter):
- The cut cells split into regions, the regions closest to the
user-defined points are retained.
Uses maxDistance for additional control.
proximity (post-filter):
- Checks the resulting faces against the original search surface
and rejects faces with a distance greater than absProximity.
ENH: restructure distance surface geometric filtering
- prefilter cells, which can be used to adjust the distance
calculation in the far field to the real distance
(not the normal distance).
This can also be used to artificially sharpen the transition
between near/far regions, if required in the future.
---
src/sampling/Make/files | 1 +
.../distanceSurface/sampledDistanceSurface.H | 16 +-
.../surface/distanceSurface/distanceSurface.C | 511 ++++++++++++++----
.../surface/distanceSurface/distanceSurface.H | 163 +++++-
.../distanceSurface/distanceSurfaceFilter.C | 421 +++++++++++++++
5 files changed, 990 insertions(+), 122 deletions(-)
create mode 100644 src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index f7d9f7d4b49..22c32f0eeb1 100644
--- a/src/sampling/Make/files
+++ b/src/sampling/Make/files
@@ -28,6 +28,7 @@ surface/cutting/cuttingSurfaceCuts.C
surface/cutting/cuttingSurfaceBase.C
surface/cutting/cuttingSurfaceBaseSelection.C
surface/distanceSurface/distanceSurface.C
+surface/distanceSurface/distanceSurfaceFilter.C
surface/isoSurface/isoSurfaceBase.C
surface/isoSurface/isoSurfaceBaseNew.C
surface/isoSurface/isoSurfaceParams.C
diff --git a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
index e3d5df2509e..ce208cb2591 100644
--- a/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/sampledDistanceSurface.H
@@ -40,6 +40,9 @@ Usage
surface1
{
type distanceSurface;
+ surfaceType triSurfaceMesh;
+ surfaceName something.obj;
+ topology proximity;
}
}
\endverbatim
@@ -48,18 +51,23 @@ Usage
\table
Property | Description | Required | Default
type | distanceSurface | yes |
- distance | Distance from surface | yes |
- signed | Use sign when distance is positive | partly |
- isoMethod | Iso-algorithm (cell/topo/point) | no | topo
- regularise | Point snapping for iso-surface | no | true
+ distance | distance from surface | no | 0
+ signed | Use sign when distance is positive | no | true
+ isoMethod | Iso-algorithm (cell/topo/point) | no | default
+ regularise | Face simplification (enum or bool) | no | true
average | Cell values from averaged point values | no | false
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
+ 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
\endtable
SourceFiles
sampledDistanceSurface.C
+ sampledDistanceSurfaceTemplates.C
\*---------------------------------------------------------------------------*/
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.C b/src/sampling/surface/distanceSurface/distanceSurface.C
index 0c6f4423ff3..be57d3b122b 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.C
+++ b/src/sampling/surface/distanceSurface/distanceSurface.C
@@ -41,6 +41,205 @@ namespace Foam
}
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::Enum
+<
+ Foam::distanceSurface::topologyFilterType
+>
+Foam::distanceSurface::topoFilterNames_
+({
+ { topologyFilterType::NONE, "none" },
+ { topologyFilterType::LARGEST_REGION, "largestRegion" },
+ { topologyFilterType::NEAREST_POINTS, "nearestPoints" },
+ { topologyFilterType::PROXIMITY, "proximity" },
+});
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Check that all point hits are valid
+static inline void checkAllHits(const UList<pointIndexHit>& nearest)
+{
+ label notHit = 0;
+ for (const pointIndexHit& pHit : nearest)
+ {
+ if (!pHit.hit())
+ {
+ ++notHit;
+ }
+ }
+
+ if (notHit)
+ {
+ FatalErrorInFunction
+ << "Had " << notHit << " from " << nearest.size()
+ << " without a point hit" << endl
+ << abort(FatalError);
+ }
+}
+
+
+// Normal distance from surface hit point to a point in the mesh
+static inline scalar normalDistance_zero
+(
+ const point& pt,
+ const pointIndexHit& pHit,
+ const vector& norm
+)
+{
+ const vector diff(pt - pHit.point());
+
+ return (diff & norm);
+}
+
+
+// Signed distance from surface hit point to a point in the mesh,
+// the sign is dictated by the normal
+static inline scalar normalDistance_nonzero
+(
+ const point& pt,
+ const pointIndexHit& pHit,
+ const vector& norm
+)
+{
+ const vector diff(pt - pHit.point());
+ const scalar normDist = (diff & norm);
+
+ return Foam::sign(normDist) * Foam::mag(diff);
+}
+
+
+// Normal distance from surface hit point to a point in the mesh
+static inline void calcNormalDistance_zero
+(
+ scalarField& distance,
+ const pointField& points,
+ const List<pointIndexHit>& nearest,
+ const vectorField& normals
+)
+{
+ forAll(nearest, i)
+ {
+ distance[i] =
+ normalDistance_zero(points[i], nearest[i], normals[i]);
+ }
+}
+
+
+// Signed distance from surface hit point -> point in the mesh,
+// the sign is dictated by the normal
+static inline void calcNormalDistance_nonzero
+(
+ scalarField& distance,
+ const pointField& points,
+ const List<pointIndexHit>& nearest,
+ const vectorField& normals
+)
+{
+ forAll(nearest, i)
+ {
+ distance[i] =
+ normalDistance_nonzero(points[i], nearest[i], normals[i]);
+ }
+}
+
+
+// Close to the surface: normal distance from surface hit point
+// Far from surface: distance from surface hit point
+//
+// Note
+// This switch may be helpful when working directly with
+// distance/gradient fields. Has low overhead otherwise.
+// May be replaced in the future (2020-11)
+static inline void calcNormalDistance_filtered
+(
+ scalarField& distance,
+ const bitSet& ignoreLocation,
+ const pointField& points,
+ const List<pointIndexHit>& nearest,
+ const vectorField& normals
+)
+{
+ forAll(nearest, i)
+ {
+ if (ignoreLocation.test(i))
+ {
+ distance[i] =
+ normalDistance_nonzero(points[i], nearest[i], normals[i]);
+ }
+ else
+ {
+ distance[i] =
+ normalDistance_zero(points[i], nearest[i], normals[i]);
+ }
+ }
+}
+
+
+// Flat surfaces (eg, a plane) have an extreme change in
+// 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
+template<bool WantPointFilter = false>
+static bitSet simpleGeometricFilter
+(
+ bitSet& ignoreCells,
+ const List<pointIndexHit>& nearest,
+ const polyMesh& mesh
+)
+{
+ // A deny filter. Initially false (accept everything)
+ ignoreCells.resize(mesh.nCells());
+
+ bitSet pointFilter;
+ if (WantPointFilter)
+ {
+ // Create as accept filter. Initially false (deny everything)
+ pointFilter.resize(mesh.nPoints());
+ }
+
+ boundBox cellBb;
+
+ forAll(nearest, celli)
+ {
+ const point& pt = nearest[celli].point();
+
+ const labelList& cPoints = mesh.cellPoints(celli);
+
+ cellBb.clear();
+ cellBb.add(mesh.points(), cPoints);
+
+ // Expand slightly to catch corners
+ cellBb.inflate(0.1);
+
+ if (!cellBb.contains(pt))
+ {
+ ignoreCells.set(celli);
+ }
+ else if (WantPointFilter)
+ {
+ // Good cell candidate, accept its points
+ pointFilter.set(cPoints);
+ }
+ }
+
+ // Flip from accept to deny filter
+ pointFilter.flip();
+
+ return pointFilter;
+}
+
+
+} // End namespace Foam
+
+
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::distanceSurface::distanceSurface
@@ -51,7 +250,7 @@ Foam::distanceSurface::distanceSurface
)
:
mesh_(mesh),
- surfPtr_
+ geometryPtr_
(
searchableSurface::New
(
@@ -68,10 +267,13 @@ Foam::distanceSurface::distanceSurface
dict
)
),
- distance_(dict.get<scalar>("distance")),
- signed_
+ distance_(dict.getOrDefault<scalar>("distance", 0)),
+ withZeroDistance_(equal(distance_, 0)),
+ withSignDistance_
(
- distance_ < 0 || equal(distance_, Zero) || dict.get<bool>("signed")
+ withZeroDistance_
+ || (distance_ < 0)
+ || dict.getOrDefault<bool>("signed", true)
),
isoParams_
(
@@ -79,9 +281,55 @@ Foam::distanceSurface::distanceSurface
isoSurfaceParams::ALGO_TOPO,
isoSurfaceParams::filterType::DIAGCELL
),
+ topoFilter_
+ (
+ topoFilterNames_.getOrDefault
+ (
+ "topology",
+ dict,
+ topologyFilterType::NONE
+ )
+ ),
+ nearestPoints_(),
+ maxDistanceSqr_(Foam::sqr(GREAT)),
+ absProximity_(dict.getOrDefault<scalar>("absProximity", 1e-5)),
+ cellDistancePtr_(nullptr),
+ pointDistance_(),
surface_(),
meshCells_(),
isoSurfacePtr_(nullptr)
+{
+ if (topologyFilterType::NEAREST_POINTS == topoFilter_)
+ {
+ dict.readEntry("nearestPoints", nearestPoints_);
+ }
+
+ if (dict.readIfPresent("maxDistance", maxDistanceSqr_))
+ {
+ maxDistanceSqr_ = Foam::sqr(maxDistanceSqr_);
+ }
+}
+
+
+Foam::distanceSurface::distanceSurface
+(
+ const polyMesh& mesh,
+ const word& surfaceType,
+ const word& surfaceName,
+ const isoSurfaceParams& params,
+ const bool interpolate
+)
+:
+ distanceSurface
+ (
+ mesh,
+ interpolate,
+ surfaceType,
+ surfaceName,
+ scalar(0),
+ true, // redundant - must be signed
+ params
+ )
{}
@@ -97,7 +345,7 @@ Foam::distanceSurface::distanceSurface
)
:
mesh_(mesh),
- surfPtr_
+ geometryPtr_
(
searchableSurface::New
(
@@ -115,11 +363,20 @@ Foam::distanceSurface::distanceSurface
)
),
distance_(distance),
- signed_
+ withZeroDistance_(equal(distance_, 0)),
+ withSignDistance_
(
- useSignedDistance || distance_ < 0 || equal(distance_, Zero)
+ withZeroDistance_
+ || (distance_ < 0)
+ || useSignedDistance
),
isoParams_(params),
+ topoFilter_(topologyFilterType::NONE),
+ nearestPoints_(),
+ maxDistanceSqr_(Foam::sqr(GREAT)),
+ absProximity_(1e-5),
+ cellDistancePtr_(nullptr),
+ pointDistance_(),
surface_(),
meshCells_(),
isoSurfacePtr_(nullptr)
@@ -140,7 +397,10 @@ void Foam::distanceSurface::createGeometry()
surface_.clear();
meshCells_.clear();
- const fvMesh& fvm = static_cast<const fvMesh&>(mesh_);
+ // Doing searches on this surface
+ const searchableSurface& geom = geometryPtr_();
+
+ const fvMesh& fvmesh = static_cast<const fvMesh&>(mesh_);
// Distance to cell centres
// ~~~~~~~~~~~~~~~~~~~~~~~~
@@ -152,132 +412,126 @@ void Foam::distanceSurface::createGeometry()
IOobject
(
"distanceSurface.cellDistance",
- fvm.time().timeName(),
- fvm.time(),
+ fvmesh.time().timeName(),
+ fvmesh.time(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
- fvm,
- dimensionedScalar(dimLength, Zero)
+ fvmesh,
+ dimensionedScalar(dimLength, GREAT)
)
);
- volScalarField& cellDistance = *cellDistancePtr_;
+ auto& cellDistance = *cellDistancePtr_;
- // For distance = 0 (and isoSurfaceCell) we apply additional filtering
+
+ // For distance = 0 we apply additional geometric filtering
// to limit the extent of open edges.
+ //
+ // Does not work with ALGO_POINT
+
+ bitSet ignoreCells, ignoreCellPoints;
- const bool isZeroDist = equal(distance_, Zero);
const bool filterCells =
(
- isZeroDist
+ withZeroDistance_
&& isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
);
- bitSet ignoreCells;
- if (filterCells)
- {
- ignoreCells.resize(fvm.nCells());
- }
// Internal field
{
- const pointField& cc = fvm.C();
+ const pointField& cc = fvmesh.C();
scalarField& fld = cellDistance.primitiveFieldRef();
List<pointIndexHit> nearest;
- surfPtr_().findNearest
+ geom.findNearest
(
cc,
- scalarField(cc.size(), GREAT),
+ // Use initialized field (GREAT) to limit search too
+ fld,
nearest
);
+ checkAllHits(nearest);
- if (signed_ || isZeroDist)
+ // Geometric pre-filtering when distance == 0
+ if (filterCells)
{
- vectorField norms;
- surfPtr_().getNormal(nearest, norms);
+ ignoreCellPoints =
+ simpleGeometricFilter<false>(ignoreCells, nearest, fvmesh);
+ }
- boundBox cellBb;
+ if (withSignDistance_)
+ {
+ vectorField norms;
+ geom.getNormal(nearest, norms);
- forAll(norms, i)
+ if (filterCells)
{
- const point diff(cc[i] - nearest[i].hitPoint());
-
- fld[i] =
+ // With inside/outside switching (see note above)
+ calcNormalDistance_filtered
(
- isZeroDist // Use normal distance
- ? (diff & norms[i])
- : Foam::sign(diff & norms[i]) * Foam::mag(diff)
+ fld,
+ ignoreCells,
+ cc,
+ nearest,
+ norms
);
-
- if (filterCells)
- {
- cellBb.clear();
- cellBb.add(fvm.points(), fvm.cellPoints(i));
-
- // Expand slightly to catch corners
- cellBb.inflate(0.1);
-
- if (!cellBb.contains(nearest[i].hitPoint()))
- {
- ignoreCells.set(i);
- }
- }
+ }
+ else if (withZeroDistance_)
+ {
+ calcNormalDistance_zero(fld, cc, nearest, norms);
+ }
+ else
+ {
+ calcNormalDistance_nonzero(fld, cc, nearest, norms);
}
}
else
{
- forAll(nearest, i)
- {
- fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
-
- // No filtering for unsigned or distance != 0.
- }
+ calcAbsoluteDistance(fld, cc, nearest);
}
}
+
// Patch fields
{
- volScalarField::Boundary& cellDistanceBf =
- cellDistance.boundaryFieldRef();
-
- forAll(fvm.C().boundaryField(), patchi)
+ forAll(fvmesh.C().boundaryField(), patchi)
{
- const pointField& cc = fvm.C().boundaryField()[patchi];
- fvPatchScalarField& fld = cellDistanceBf[patchi];
+ const pointField& cc = fvmesh.C().boundaryField()[patchi];
+ scalarField& fld = cellDistance.boundaryFieldRef()[patchi];
List<pointIndexHit> nearest;
- surfPtr_().findNearest
+ geom.findNearest
(
cc,
scalarField(cc.size(), GREAT),
nearest
);
+ checkAllHits(nearest);
- if (signed_)
+ if (withSignDistance_)
{
vectorField norms;
- surfPtr_().getNormal(nearest, norms);
+ geom.getNormal(nearest, norms);
+
+ if (withZeroDistance_)
+ {
+ // Slight inconsistency in boundary vs interior when
+ // cells are filtered, but the patch fields are only
+ // used by isoSurfacePoint, and filtering is disabled
+ // for that anyhow.
- forAll(norms, i)
+ calcNormalDistance_zero(fld, cc, nearest, norms);
+ }
+ else
{
- const point diff(cc[i] - nearest[i].hitPoint());
-
- fld[i] =
- (
- isZeroDist // Use normal distance
- ? (diff & norms[i])
- : Foam::sign(diff & norms[i]) * Foam::mag(diff)
- );
+ calcNormalDistance_nonzero(fld, cc, nearest, norms);
}
}
else
{
- forAll(nearest, i)
- {
- fld[i] = Foam::mag(cc[i] - nearest[i].hitPoint());
- }
+ calcAbsoluteDistance(fld, cc, nearest);
}
}
}
@@ -288,44 +542,86 @@ void Foam::distanceSurface::createGeometry()
// Distance to points
- pointDistance_.resize(fvm.nPoints());
+ pointDistance_.resize(fvmesh.nPoints());
+ pointDistance_ = GREAT;
{
- const pointField& pts = fvm.points();
+ const pointField& pts = fvmesh.points();
+ scalarField& fld = pointDistance_;
List<pointIndexHit> nearest;
- surfPtr_().findNearest
+ geom.findNearest
(
pts,
- scalarField(pts.size(), GREAT),
+ // Use initialized field (GREAT) to limit search too
+ pointDistance_,
nearest
);
+ checkAllHits(nearest);
- if (signed_)
+ if (withSignDistance_)
{
vectorField norms;
- surfPtr_().getNormal(nearest, norms);
+ geom.getNormal(nearest, norms);
- forAll(norms, i)
+ if (filterCells)
{
- const point diff(pts[i] - nearest[i].hitPoint());
-
- pointDistance_[i] =
+ // With inside/outside switching (see note above)
+ calcNormalDistance_filtered
(
- isZeroDist // Use normal distance
- ? (diff & norms[i])
- : Foam::sign(diff & norms[i]) * Foam::mag(diff)
+ fld,
+ ignoreCellPoints,
+ pts,
+ nearest,
+ norms
);
}
+ else if (withZeroDistance_)
+ {
+ calcNormalDistance_zero(fld, pts, nearest, norms);
+ }
+ else
+ {
+ calcNormalDistance_nonzero(fld, pts, nearest, norms);
+ }
}
else
{
- forAll(nearest, i)
- {
- pointDistance_[i] = Foam::mag(pts[i] - nearest[i].hitPoint());
- }
+ calcAbsoluteDistance(fld, pts, nearest);
+ }
+ }
+
+
+ // Don't need ignoreCells if there is nothing to ignore.
+ if (ignoreCells.none())
+ {
+ ignoreCells.clearStorage();
+ }
+ else if (filterCells && topologyFilterType::NONE != topoFilter_)
+ {
+ // For refine blocked cells (eg, checking actual cells cut)
+ isoSurfaceBase isoCutter
+ (
+ mesh_,
+ cellDistance,
+ pointDistance_,
+ distance_
+ );
+
+ if (topologyFilterType::LARGEST_REGION == topoFilter_)
+ {
+ refineBlockedCells(ignoreCells, isoCutter);
+ filterKeepLargestRegion(ignoreCells);
+ }
+ else if (topologyFilterType::NEAREST_POINTS == topoFilter_)
+ {
+ refineBlockedCells(ignoreCells, isoCutter);
+ filterKeepNearestRegions(ignoreCells);
}
}
+ // Don't need point filter beyond this point
+ ignoreCellPoints.clearStorage();
+
if (debug)
{
@@ -336,13 +632,13 @@ void Foam::distanceSurface::createGeometry()
IOobject
(
"distanceSurface.pointDistance",
- fvm.time().timeName(),
- fvm.time(),
+ fvmesh.time().timeName(),
+ fvmesh.time(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
- pointMesh::New(fvm),
+ pointMesh::New(fvmesh),
dimensionedScalar(dimLength, Zero)
);
pDist.primitiveFieldRef() = pointDistance_;
@@ -351,13 +647,6 @@ void Foam::distanceSurface::createGeometry()
pDist.write();
}
- // Don't need ignoreCells if there is nothing to ignore.
- if (!ignoreCells.any())
- {
- ignoreCells.clear();
- }
-
-
isoSurfacePtr_.reset
(
isoSurfaceBase::New
@@ -370,9 +659,16 @@ void Foam::distanceSurface::createGeometry()
)
);
+
// ALGO_POINT still needs cell, point fields (for interpolate)
// The others can do straight transfer
- if (isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT)
+
+ // But also flatten into a straight transfer for proximity filtering
+ if
+ (
+ isoParams_.algorithm() != isoSurfaceParams::ALGO_POINT
+ || topologyFilterType::PROXIMITY == topoFilter_
+ )
{
surface_.transfer(static_cast<meshedSurface&>(*isoSurfacePtr_));
meshCells_.transfer(isoSurfacePtr_->meshCells());
@@ -382,6 +678,11 @@ void Foam::distanceSurface::createGeometry()
pointDistance_.clear();
}
+ if (topologyFilterType::PROXIMITY == topoFilter_)
+ {
+ filterByProximity();
+ }
+
if (debug)
{
print(Pout);
diff --git a/src/sampling/surface/distanceSurface/distanceSurface.H b/src/sampling/surface/distanceSurface/distanceSurface.H
index b56f49375b3..22738393291 100644
--- a/src/sampling/surface/distanceSurface/distanceSurface.H
+++ b/src/sampling/surface/distanceSurface/distanceSurface.H
@@ -32,27 +32,82 @@ Description
Uses an iso-surface algorithm (cell, topo, point) for constructing the
distance surface.
+ For a zero-distance surface, it performs additional checks and supports
+ filtering to handle the surface boundaries.
+
Usage
+ Example of function object partial specification:
+ \verbatim
+ surfaces
+ {
+ surface1
+ {
+ type distanceSurface;
+ surfaceType triSurfaceMesh;
+ surfaceName something.obj;
+ topology proximity;
+ }
+
+ surface2
+ {
+ type distanceSurface;
+ surfaceType triSurfaceMesh;
+ surfaceName other.obj;
+
+ topology nearestPoints;
+ nearestPoints
+ (
+ (0 0 0)
+ (10 10 0)
+ );
+
+ // Max search distance for nearestPoints
+ maxDistance 0.005;
+ }
+ }
+ \endverbatim
+
Dictionary controls:
\table
Property | Description | Required | Default
- distance | distance from surface | yes |
- signed | Use sign when distance is positive | partly |
- isoMethod | Iso-algorithm (cell/topo/point) | no | topo
+ distance | distance from surface | no | 0
+ signed | Use sign when distance is positive | no | true
+ isoMethod | Iso-algorithm (cell/topo/point) | no | default
regularise | Face simplification (enum or bool) | no | true
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
+ 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)
+
+ - \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
For distance = 0, some special adjustments.
- Always signed (ignoring the input value).
- Use normal distance from surface (for better treatment of open edges).
- - When 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.
+ - Additional checks for open surfaces edges are used to limit the extend
+ of resulting distance surface.
+ The resulting surface elements will, however, contain partial cell
+ coverage. NB: Not applicable if the \c point isoMethod is used.
The keyword \c cell (bool value) which was use in 1906 and earlier to switch
between point/cell algorithms is now ignored (2020-12).
@@ -70,6 +125,7 @@ SourceFiles
#include "sampledSurface.H"
#include "searchableSurface.H"
#include "isoSurfaceBase.H"
+#include "pointIndexHit.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -82,23 +138,53 @@ namespace Foam
class distanceSurface
{
+ // Data Types
+
+ //- The type of pre/post face-filtering
+ enum class topologyFilterType : uint8_t
+ {
+ NONE, //!< No additional filtering
+ LARGEST_REGION, //!< Retain largest region
+ NEAREST_POINTS, //!< Retain regions nearest to the points
+ PROXIMITY //!< Post-filter for surface proximity
+ };
+
+ //- Names for the topology filter
+ static const Enum<topologyFilterType> topoFilterNames_;
+
+
// Private Data
//- Reference to mesh
const polyMesh& mesh_;
- //- Surface
- const autoPtr<searchableSurface> surfPtr_;
+ //- Searchable surface
+ const autoPtr<searchableSurface> geometryPtr_;
//- Distance value
const scalar distance_;
- //- Signed distance
- const bool signed_;
+ //- Distance is zero. Implies signed and additional optimizations
+ const bool withZeroDistance_;
+
+ //- Use signed distance
+ const bool withSignDistance_;
//- Parameters for iso-surface (algorithm, filter, mergeTol, etc)
isoSurfaceParams isoParams_;
+ //- Optional topology face-filtering
+ topologyFilterType topoFilter_;
+
+ //- Points for nearest-points segmentation
+ pointField nearestPoints_;
+
+ //- Max search distance squared (for nearestPoints)
+ scalar maxDistanceSqr_;
+
+ //- Max distance for proximity check (post-filtering)
+ scalar absProximity_;
+
//- Distance to cell centres
autoPtr<volScalarField> cellDistancePtr_;
@@ -118,6 +204,24 @@ class distanceSurface
mutable autoPtr<isoSurfaceBase> isoSurfacePtr_;
+ // Private Member Functions
+
+ //- Absolute distances from hit points
+ // Hit/miss checks have been done elsewhere.
+ static inline void calcAbsoluteDistance
+ (
+ scalarField& distance,
+ const pointField& points,
+ const List<pointIndexHit>& nearest
+ )
+ {
+ forAll(nearest, i)
+ {
+ distance[i] = Foam::mag(points[i] - nearest[i].point());
+ }
+ }
+
+
protected:
// Protected Member Functions
@@ -145,6 +249,29 @@ protected:
}
+ // Private Member Functions
+
+ //- Re-filter the blocked cells based on the anticipated cuts
+ // Uses a lightweight variant of cutting.
+ bool refineBlockedCells
+ (
+ bitSet& ignoreCells,
+ const isoSurfaceBase& isoContext
+ ) const;
+
+ //- 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;
+
+
public:
//- Runtime type information
@@ -161,6 +288,16 @@ public:
const dictionary& dict
);
+ //- Construct from components with zero-distanced
+ distanceSurface
+ (
+ const polyMesh& mesh,
+ const word& surfaceType,
+ const word& surfaceName,
+ const isoSurfaceParams& params = isoSurfaceParams(),
+ const bool interpolate = false
+ );
+
//- Construct from components
distanceSurface
(
@@ -186,11 +323,11 @@ public:
//- The name of the underlying searchableSurface
const word& surfaceName() const
{
- return surfPtr_->name();
+ return geometryPtr_->name();
}
//- The distance to the underlying searchableSurface
- scalar distance() const
+ scalar distance() const noexcept
{
return distance_;
}
diff --git a/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C b/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
new file mode 100644
index 00000000000..95990583024
--- /dev/null
+++ b/src/sampling/surface/distanceSurface/distanceSurfaceFilter.C
@@ -0,0 +1,421 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "distanceSurface.H"
+#include "regionSplit.H"
+#include "syncTools.H"
+#include "vtkSurfaceWriter.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+bool Foam::distanceSurface::refineBlockedCells
+(
+ bitSet& ignoreCells,
+ const isoSurfaceBase& isoCutter
+) const
+{
+ // With the cell/point distance fields we can take a second pass at
+ // pre-filtering.
+ // This duplicates how cut detection is determined in the cell/topo
+ // algorithms but is fairly inexpensive (creates no geometry)
+
+ bool changed = false;
+
+ for (label celli = 0; celli < mesh_.nCells(); ++celli)
+ {
+ if (ignoreCells.test(celli))
+ {
+ continue;
+ }
+
+ auto cut = isoCutter.getCellCutType(celli);
+ if (!(cut & isoSurfaceBase::ANYCUT))
+ {
+ ignoreCells.set(celli);
+ changed = true;
+ }
+ }
+
+ return changed;
+}
+
+
+Foam::bitSet Foam::distanceSurface::filterPrepareRegionSplit
+(
+ const bitSet& ignoreCells
+) const
+{
+ // Prepare for region split
+
+ bitSet blockedFaces(mesh_.nFaces());
+
+ const labelList& faceOwn = mesh_.faceOwner();
+ const labelList& faceNei = mesh_.faceNeighbour();
+
+ // Could be more efficient
+ for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
+ {
+ // If only one cell is blocked, the face corresponds
+ // to an exposed subMesh face
+
+ if
+ (
+ ignoreCells.test(faceOwn[facei])
+ != ignoreCells.test(faceNei[facei])
+ )
+ {
+ blockedFaces.set(facei);
+ }
+ }
+
+ for (const polyPatch& patch : mesh_.boundaryMesh())
+ {
+ if (!patch.coupled())
+ {
+ continue;
+ }
+
+ forAll(patch, patchFacei)
+ {
+ const label facei = patch.start() + patchFacei;
+ if (ignoreCells.test(faceOwn[facei]))
+ {
+ blockedFaces.set(facei);
+ }
+ }
+ }
+
+ syncTools::syncFaceList(mesh_, blockedFaces, xorEqOp<unsigned int>());
+
+ return blockedFaces;
+}
+
+
+void Foam::distanceSurface::filterKeepLargestRegion
+(
+ bitSet& ignoreCells
+) const
+{
+ // For region split
+ bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells));
+
+ // Split region
+ regionSplit rs(mesh_, blockedFaces);
+ blockedFaces.clearStorage();
+
+ const labelList& regionColour = rs;
+
+ // Identical number of regions on all processors
+ labelList nCutsPerRegion(rs.nRegions(), Zero);
+
+ // Count cut cells (ie, unblocked)
+ forAll(regionColour, celli)
+ {
+ if (!ignoreCells.test(celli))
+ {
+ ++nCutsPerRegion[regionColour[celli]];
+ }
+ }
+
+ // Sum totals from all processors
+ Pstream::listCombineGather(nCutsPerRegion, plusEqOp<label>());
+
+
+ // Define which regions to keep
+ boolList keepRegion(rs.nRegions(), false);
+
+ if (Pstream::master())
+ {
+ const label largest = findMax(nCutsPerRegion);
+
+ if (largest == -1)
+ {
+ // Should not happen
+ keepRegion = true;
+ }
+ else
+ {
+ keepRegion[largest] = true;
+ }
+
+ if (debug)
+ {
+ Info<< "Had " << sum(nCutsPerRegion) << " cuts, in "
+ << nCutsPerRegion.size() << " regions, largest is "
+ << largest << ": " << flatOutput(nCutsPerRegion) << nl;
+ }
+ }
+
+ Pstream::scatter(keepRegion);
+
+ forAll(regionColour, celli)
+ {
+ if (!keepRegion.test(regionColour[celli]))
+ {
+ ignoreCells.set(celli);
+ }
+ }
+}
+
+
+void Foam::distanceSurface::filterKeepNearestRegions
+(
+ bitSet& ignoreCells
+) const
+{
+ if (nearestPoints_.empty())
+ {
+ WarningInFunction
+ << "Ignoring nearestPoints - no points provided" << nl
+ << endl;
+ return;
+ }
+
+ // For region split
+ bitSet blockedFaces(filterPrepareRegionSplit(ignoreCells));
+
+ // Split region
+ regionSplit rs(mesh_, blockedFaces);
+ blockedFaces.clearStorage();
+
+ const labelList& regionColour = rs;
+
+ const pointField& cc = mesh_.cellCentres();
+ const pointField& nearPts = nearestPoints_;
+
+ // The magSqr distance and region
+ typedef Tuple2<scalar, label> nearInfo;
+ List<nearInfo> nearest(nearPts.size(), nearInfo(GREAT, -1));
+
+ // Also collect cuts per region, may be useful for rejecting
+ // small regions. Code as per filterKeepLargestRegion
+ labelList nCutsPerRegion(rs.nRegions(), Zero);
+
+ forAll(cc, celli)
+ {
+ if (ignoreCells.test(celli))
+ {
+ continue;
+ }
+
+ const point& pt = cc[celli];
+ const label regioni = regionColour[celli];
+
+ ++nCutsPerRegion[regioni];
+
+ label pointi = 0;
+ for (nearInfo& near : nearest)
+ {
+ const scalar distSqr = magSqr(nearPts[pointi] - pt);
+ ++pointi;
+
+ if (distSqr < near.first())
+ {
+ near.first() = distSqr;
+ near.second() = regioni;
+ }
+ }
+ }
+
+ // Sum totals from all processors
+ Pstream::listCombineGather(nCutsPerRegion, plusEqOp<label>());
+
+ // Get nearest
+ Pstream::listCombineGather(nearest, minFirstEqOp<scalar>());
+
+
+ // Define which regions to keep
+
+ boolList keepRegion(rs.nRegions(), false);
+
+ if (Pstream::master())
+ {
+ const label largest = findMax(nCutsPerRegion);
+
+ for (const nearInfo& near : nearest)
+ {
+ const scalar distSqr = near.first();
+ const label regioni = near.second();
+
+ if (regioni != -1 && distSqr < maxDistanceSqr_)
+ {
+ keepRegion[regioni] = true;
+ }
+ }
+
+ if (debug)
+ {
+ Info<< "Had " << sum(nCutsPerRegion) << " cuts, in "
+ << nCutsPerRegion.size() << " regions, largest is "
+ << largest << ": " << flatOutput(nCutsPerRegion) << nl;
+
+ Info<< "nearestPoints (max distance = "
+ << sqrt(maxDistanceSqr_) << ")" << nl;
+
+ forAll(nearPts, pointi)
+ {
+ const scalar dist = sqrt(nearest[pointi].first());
+ const label regioni = nearest[pointi].second();
+
+ Info<< " " << nearPts[pointi] << " region "
+ << regioni << " distance "
+ << dist;
+
+ if (!keepRegion.test(regioni))
+ {
+ Info<< " too far";
+ }
+ Info<< nl;
+ }
+ }
+ }
+
+ Pstream::scatter(keepRegion);
+
+ forAll(regionColour, celli)
+ {
+ if (!keepRegion.test(regionColour[celli]))
+ {
+ ignoreCells.set(celli);
+ }
+ }
+}
+
+
+void Foam::distanceSurface::filterByProximity()
+{
+ const searchableSurface& geom = geometryPtr_();
+
+ // Filtering for faces
+ const pointField& fc = surface_.faceCentres();
+
+ bitSet faceSelection(fc.size());
+ label nTrimmed = 0;
+
+
+ // For each face
+ scalarField faceDistance(fc.size(), GREAT);
+ scalarField faceNormalDistance; // Debugging only
+ {
+ List<pointIndexHit> nearest;
+ geom.findNearest
+ (
+ fc,
+ // Use initialized field (GREAT) to limit search too
+ faceDistance,
+ nearest
+ );
+ calcAbsoluteDistance(faceDistance, fc, nearest);
+
+ // Heavier debugging
+ if (debug & 4)
+ {
+ vectorField norms;
+ geom.getNormal(nearest, norms);
+
+ faceNormalDistance.resize(fc.size());
+
+ forAll(nearest, i)
+ {
+ const vector diff(fc[i] - nearest[i].point());
+
+ faceNormalDistance[i] = Foam::mag((diff & norms[i]));
+ }
+ }
+ }
+
+ // Note
+ // Proximity checks using the face points (nearest hit) to establish
+ // a length scale are too fragile. Can easily have stretched faces
+ // where the centre is less than say 0.3-0.5 of the centre-point distance
+ // but they are still outside.
+
+ // Using the absolute proximity of the face centres is more robust.
+
+
+ // Consider the absolute proximity of the face centres
+ forAll(faceDistance, facei)
+ {
+ if (faceDistance[facei] <= absProximity_)
+ {
+ faceSelection.set(facei);
+ }
+ else
+ {
+ ++nTrimmed;
+
+ if (debug & 2)
+ {
+ Pout<< "trim reject: "
+ << faceDistance[facei] << nl;
+ }
+ }
+ }
+
+
+ // Heavier debugging
+ if (debug & 4)
+ {
+ labelField faceFilterStatus(faceSelection.size(), Zero);
+
+ for (const label facei : faceSelection)
+ {
+ faceFilterStatus[facei] = 1;
+ }
+
+ const fileName outputName(surfaceName() + "-proximity-filter");
+
+ Info<< "Writing debug surface: " << outputName << nl;
+
+ surfaceWriters::vtkWriter writer
+ (
+ surface_.points(),
+ surface_, // faces
+ outputName
+ );
+
+ writer.write("absolute-distance", faceDistance);
+ writer.write("normal-distance", faceNormalDistance);
+ writer.write("filter-state", faceFilterStatus);
+ }
+
+
+ if (returnReduce(nTrimmed, sumOp<label>()) != 0)
+ {
+ labelList pointMap, faceMap;
+ meshedSurface filtered
+ (
+ surface_.subsetMesh(faceSelection, pointMap, faceMap)
+ );
+ surface_.transfer(filtered);
+
+ meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
+ }
+}
+
+
+// ************************************************************************* //
--
GitLab