diff --git a/src/sampling/Make/files b/src/sampling/Make/files
index f7d9f7d4b490c37da7dc6ae0e3dd02ba12c13d86..22c32f0eeb11f19d3686dc45e44ecb79d113178e 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 e3d5df2509e4e0eb9294a4fc79f72c9ee756801a..ce208cb2591ea72c4c7944c8def26bbd69949ee7 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 0c6f4423ff37f7740da8c04f8b689fe19f5069d7..be57d3b122b6aa294f6f73a6bc65c580d2f07265 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 b56f49375b3362302bde2eb44bb90aff8771ded7..22738393291b3783e2b36059adca04144f461a9c 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 0000000000000000000000000000000000000000..9599058302400ac6dc8ee707ca9b6996880841fd
--- /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)();
+ }
+}
+
+
+// ************************************************************************* //