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

ENH: additional filtering for distance surface (#1950)

Browse Source 
- 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.
This commit is contained in:
Mark Olesen
2020-12-09 11:36:49 +01:00
committed by Andrew Heather
parent 7a99866db0
commit 4c3a95c808
5 changed files with 989 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/sampling/Make/files
View File

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

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

@ -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
\*---------------------------------------------------------------------------*/

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

@ -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)
{
ignoreCellPoints =
simpleGeometricFilter<false>(ignoreCells, nearest, fvmesh);
}
if (withSignDistance_)
{
vectorField norms;
surfPtr_().getNormal(nearest, norms);
geom.getNormal(nearest, norms);
boundBox cellBb;
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);
forAll(norms, i)
if (withZeroDistance_)
{
const point diff(cc[i] - nearest[i].hitPoint());
// 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.
fld[i] =
(
isZeroDist // Use normal distance
? (diff & norms[i])
: Foam::sign(diff & norms[i]) * Foam::mag(diff)
);
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());
}
calcAbsoluteDistance(fld, cc, nearest);
}
}
}
@ -288,45 +542,87 @@ 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)
{
Pout<< "Writing cell distance:" << cellDistance.objectPath() << endl;
@ -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);

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

@ -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_;
}

421
src/sampling/surface/distanceSurface/distanceSurfaceFilter.C Normal file
View File

@ -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)();
}
}
// ************************************************************************* //