From 4fa63d9bc6c6ec7c74b1c101e5454ce02a1e1de0 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Fri, 19 Jul 2019 13:33:26 +0200
Subject: [PATCH] ENH: improvements for isoSurfaceTopo erosion (#1374)
- adapted openfoam.org code. Original commit message:
Instead of adapting tet base points cell-by-cell, the dangling
points are pre-computed and then the adaptations to the base points
are made face-by-face. This correctly adapts faces which have
different dangling points relative to the owner and neighbour cells.
---
.../surface/isoSurface/isoSurfaceTopo.C | 164 +++++++++---------
1 file changed, 84 insertions(+), 80 deletions(-)
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
index 744d995c7b..f2fddbaa1b 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -254,126 +254,130 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
// Pre-filter: mark all cells with illegal base points
- labelHashSet problemCells(cells.size()/128);
+ bitSet problemCells(cells.size());
+
forAll(tetBasePtIs_, facei)
{
if (tetBasePtIs_[facei] == -1)
{
- problemCells.insert(faceOwner[facei]);
+ problemCells.set(faceOwner[facei]);
+
if (mesh_.isInternalFace(facei))
{
- problemCells.insert(faceNeighbour[facei]);
+ problemCells.set(faceNeighbour[facei]);
}
}
}
- label nAdapted = 0;
-
+ // Mark all points that are shared by just two faces within an adjacent
+ // problem cell as problematic
+ bitSet problemPoints(mesh_.points().size());
- // Number of times a point occurs in a cell. Used to detect dangling
- // vertices (count = 2)
- Map<label> pointCount;
-
- // Analyse problem cells for points shared by two faces only
- forAll(cells, celli)
{
- if (problemCells.found(celli))
- {
- const cell& cFaces = cells[celli];
+ // Number of times a point occurs in a cell.
+ // Used to detect dangling vertices (count = 2)
+ Map<label> pointCount;
+ // Analyse problem cells for points shared by two faces only
+ for (const label celli : problemCells)
+ {
pointCount.clear();
- forAll(cFaces, i)
+ for (const label facei : cells[celli])
{
- const label facei = cFaces[i];
- const face& f = faces[facei];
- forAll(f, fp)
+ for (const label pointi : faces[facei])
{
- const label pointi = f[fp];
-
- Map<label>::iterator pointFnd = pointCount.find(pointi);
- if (pointFnd == pointCount.end())
- {
- pointCount.insert(pointi, 1);
- }
- else
- {
- ++pointFnd();
- }
+ ++pointCount(pointi);
}
}
- // Check for any points with count 2
- bool haveDangling = false;
forAllConstIters(pointCount, iter)
{
- if (iter() == 1)
+ if (iter.val() == 1)
{
- FatalErrorInFunction << "point:" << iter.key()
+ FatalErrorInFunction
+ << "point:" << iter.key()
<< " at:" << mesh_.points()[iter.key()]
- << " only used by one face" << exit(FatalError);
+ << " only used by one face" << nl
+ << exit(FatalError);
}
- else if (iter() == 2)
+ else if (iter.val() == 2)
{
- haveDangling = true;
- break;
+ problemPoints.set(iter.key());
}
}
+ }
+ }
- if (haveDangling)
- {
- // Any point next to a dangling point should not be used
- // as the fan base since this would cause two duplicate
- // triangles.
- forAll(cFaces, i)
- {
- const label facei = cFaces[i];
- if (tetBasePtIs_[facei] == -1)
- {
- const face& f = faces[facei];
- // All the possible base points cause negative tets.
- // Choose the least-worst one
- scalar maxQ = -GREAT;
- label maxFp = -1;
+ // For all faces which form a part of a problem-cell, check if the base
+ // point is adjacent to any problem points. If it is, re-calculate the base
+ // point so that it is not.
+ label nAdapted = 0;
+ forAll(tetBasePtIs_, facei)
+ {
+ if
+ (
+ problemCells[faceOwner[facei]]
+ || (mesh_.isInternalFace(facei) && problemCells[faceNeighbour[facei]])
+ )
+ {
+ const face& f = faces[facei];
- label prevCount = pointCount[f.last()];
- forAll(f, fp)
- {
- label nextCount = pointCount[f[f.fcIndex(fp)]];
-
- if (prevCount > 2 && nextCount > 2)
- {
- const scalar q = minTetQ(facei, fp);
- if (q > maxQ)
- {
- maxQ = q;
- maxFp = fp;
- }
- }
- prevCount = pointCount[f[fp]];
- }
+ // Check if either of the points adjacent to the base point is a
+ // problem point. If not, the existing base point can be retained.
+ const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
- if (maxFp != -1)
- {
- // Least worst base point
- tetBasePtIs_[facei] = maxFp;
- }
- else
- {
- // No point found on face that would not result
- // in some duplicate triangle. Very rare. Do what?
- tetBasePtIs_[facei] = 0;
- }
+ if
+ (
+ !problemPoints[f[f.rcIndex(fp0)]]
+ && !problemPoints[f[f.fcIndex(fp0)]]
+ )
+ {
+ continue;
+ }
- nAdapted++;
+ // A new base point is required. Pick the point that results in the
+ // least-worst tet and which is not adjacent to any problem points.
+ scalar maxQ = -GREAT;
+ label maxFp = -1;
+ forAll(f, fp)
+ {
+ if
+ (
+ !problemPoints[f[f.rcIndex(fp)]]
+ && !problemPoints[f[f.fcIndex(fp)]]
+ )
+ {
+ const scalar q = minTetQ(facei, fp);
+ if (q > maxQ)
+ {
+ maxQ = q;
+ maxFp = fp;
}
}
}
+
+ if (maxFp != -1)
+ {
+ // Success! Set the new base point
+ tetBasePtIs_[facei] = maxFp;
+ }
+ else
+ {
+ // No point was found on face that would not result in some
+ // duplicate triangle. Do what? Continue and hope? Spit an
+ // error? Silently or noisily reduce the filtering level?
+
+ tetBasePtIs_[facei] = 0;
+ }
+
+ ++nAdapted;
}
}
+
if (debug)
{
Pout<< "isoSurface : adapted starting point of triangulation on "
--
GitLab