From 61dc4839f8b52dee6a0e32dcf37135d73eb62d09 Mon Sep 17 00:00:00 2001
From: Mark Olesen <Mark.Olesen@esi-group.com>
Date: Wed, 15 Sep 2021 21:39:38 +0200
Subject: [PATCH] ENH: add point snapping to iso-surface topo algorithm (#2210)
- helps avoid the creation of small face cuts (near corners, edges)
that result in zero-size faces on output.
---
.../surface/isoSurface/isoSurfaceBase.H | 22 +-
.../surface/isoSurface/isoSurfaceParams.C | 4 +-
.../surface/isoSurface/isoSurfaceParams.H | 18 +-
.../surface/isoSurface/isoSurfaceTopo.C | 663 +++++++++++-------
.../surface/isoSurface/isoSurfaceTopo.H | 81 ++-
.../isoSurface/isoSurfaceTopoTemplates.C | 37 +-
6 files changed, 493 insertions(+), 332 deletions(-)
diff --git a/src/sampling/surface/isoSurface/isoSurfaceBase.H b/src/sampling/surface/isoSurface/isoSurfaceBase.H
index 3e29657bd11..a856f06f263 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceBase.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceBase.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -127,26 +127,6 @@ protected:
// Protected Member Functions
- //- Check for tet values above/below given (iso) value
- // Result encoded as a single integer
- inline static constexpr int getTetCutIndex
- (
- const scalar a,
- const scalar b,
- const scalar c,
- const scalar d,
- const scalar isoval
- ) noexcept
- {
- return
- (
- (a < isoval ? 1 : 0)
- | (b < isoval ? 2 : 0)
- | (c < isoval ? 4 : 0)
- | (d < isoval ? 8 : 0)
- );
- }
-
//- Count the number of cuts matching the mask type
// Checks as bitmask or as zero.
static label countCutType
diff --git a/src/sampling/surface/isoSurface/isoSurfaceParams.C b/src/sampling/surface/isoSurface/isoSurfaceParams.C
index cf9494ed5d3..5687cb3a8b4 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceParams.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceParams.C
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -137,6 +137,7 @@ Foam::isoSurfaceParams::isoSurfaceParams
:
algo_(algo),
filter_(filter),
+ snap_(true),
mergeTol_(1e-6),
clipBounds_(boundBox::invertedBox)
{}
@@ -152,6 +153,7 @@ Foam::isoSurfaceParams::isoSurfaceParams
{
algo_ = getAlgorithmType(dict, algo_);
filter_ = getFilterType(dict, filter_);
+ snap_ = dict.getOrDefault("snap", true);
dict.readIfPresent("mergeTol", mergeTol_);
dict.readIfPresent("bounds", clipBounds_);
}
diff --git a/src/sampling/surface/isoSurface/isoSurfaceParams.H b/src/sampling/surface/isoSurface/isoSurfaceParams.H
index 1b794f70595..7e0ce92002e 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceParams.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceParams.H
@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -35,6 +35,7 @@ Description
isoMethod | Algorithm (cell/topo/point/default) | no | default
regularise | Face simplification (enum or bool) | no | true
mergeTol | Point merge tolerance (cell/point) | no | 1e-6
+ snap | Point snapping (topo) | no | true
bounds | Optional clip bounds | no | inverted
\endtable
@@ -100,6 +101,9 @@ private:
//- Filtering for iso-surface faces/points
filterType filter_;
+ //- Point snapping enabled
+ bool snap_;
+
//- Merge tolerance for cell/point (default: 1e-6)
scalar mergeTol_;
@@ -186,6 +190,18 @@ public:
filter_ = fltr;
}
+ //- Get point snapping flag
+ bool snap() const noexcept
+ {
+ return snap_;
+ }
+
+ //- Set point snapping flag
+ void snap(bool on) noexcept
+ {
+ snap_ = on;
+ }
+
//- Get current merge tolerance
scalar mergeTol() const noexcept
{
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
index 515ae48f992..63a0845fff1 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -51,6 +51,112 @@ namespace Foam
}
+// Get/set snapIndex (0, 1 or 2) at given position
+
+#undef SNAP_INDEX_ENCODE
+#undef SNAP_INDEX_DECODE
+#define SNAP_INDEX_ENCODE(pos, val) (((val) << (4 + 2 * pos)))
+#define SNAP_INDEX_DECODE(pos, val) (((val) >> (4 + 2 * pos)) & 0x3)
+
+
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Check for tet values above/below given (iso) value
+// Result encoded as an integer, with possible snapping information too
+inline static int getTetCutIndex
+(
+ scalar p0,
+ scalar p1,
+ scalar p2,
+ scalar p3,
+ const scalar val,
+ const bool doSnap
+) noexcept
+{
+ int cutIndex
+ (
+ (p0 < val ? 1 : 0) // point 0
+ | (p1 < val ? 2 : 0) // point 1
+ | (p2 < val ? 4 : 0) // point 2
+ | (p3 < val ? 8 : 0) // point 3
+ );
+
+ if (doSnap && cutIndex && cutIndex != 0xF)
+ {
+ // Not all below or all
+ #undef SNAP_INDEX_CALC
+ #undef SNAP_INDEX_VALUE
+
+ // Calculate distances (for snapping)
+ p0 -= val; if (cutIndex & 1) p0 *= -1;
+ p1 -= val; if (cutIndex & 2) p1 *= -1;
+ p2 -= val; if (cutIndex & 4) p2 *= -1;
+ p3 -= val; if (cutIndex & 8) p3 *= -1;
+
+ // Snap to end if less than approx 1% of the distance.
+ // - only valid if there is also a corresponding sign change
+ #define SNAP_INDEX_VALUE(d1, d2) \
+ ((d1 * 100 < d2) ? 1 : (d2 * 100 < d1) ? 2 : 0)
+
+ // Add snap index into regular edge cut index
+ #define SNAP_INDEX_CALC(pos, d1, d2, idx1, idx2) \
+ switch (cutIndex & (idx1 | idx2)) \
+ { \
+ case idx1 : /* first below, second above */ \
+ case idx2 : /* first above, second below */ \
+ cutIndex |= SNAP_INDEX_ENCODE(pos, SNAP_INDEX_VALUE(d1, d2)); \
+ break; \
+ }
+
+ SNAP_INDEX_CALC(0, p0, p1, 1, 2); // Edge 0: 0 -> 1
+ SNAP_INDEX_CALC(1, p0, p2, 1, 4); // Edge 1: 0 -> 2
+ SNAP_INDEX_CALC(2, p0, p3, 1, 8); // Edge 2: 0 -> 3
+ SNAP_INDEX_CALC(3, p3, p1, 8, 2); // Edge 3: 3 -> 1
+ SNAP_INDEX_CALC(4, p1, p2, 2, 4); // Edge 4: 1 -> 2
+ SNAP_INDEX_CALC(5, p3, p2, 8, 4); // Edge 5: 3 -> 2
+
+ #undef SNAP_INDEX_CALC
+ #undef SNAP_INDEX_VALUE
+ }
+
+ return cutIndex;
+}
+
+
+// Append three labels to list.
+// Filter out degenerate (eg, snapped) tris. Flip face as requested
+inline static void appendTriLabels
+(
+ DynamicList<label>& verts,
+ const label a,
+ const label b,
+ const label c,
+ const bool flip // Flip normals
+)
+{
+ if (a != b && b != c && c != a)
+ {
+ verts.append(a);
+ if (flip)
+ {
+ verts.append(c);
+ verts.append(b);
+ }
+ else
+ {
+ verts.append(b);
+ verts.append(c);
+ }
+ }
+}
+
+
+} // End namespace Foam
+
+
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::scalar Foam::isoSurfaceTopo::minTetQ
@@ -239,28 +345,68 @@ void Foam::isoSurfaceTopo::fixTetBasePtIs()
}
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
Foam::label Foam::isoSurfaceTopo::generatePoint
(
- const label facei,
- const bool edgeIsDiag,
+ label facei,
+ const int pointEdgeType,
const edge& vertices,
+ // Output
DynamicList<edge>& pointToVerts,
DynamicList<label>& pointToFace,
DynamicList<bool>& pointFromDiag,
- EdgeMap<label>& vertsToPoint
-) const
+
+ // Bookkeeping
+ EdgeMap<label>& vertsToPoint,
+ Map<label>& snapVertsLookup
+)
{
// Generate new point, unless it already exists for edge
+ // or corresponds to a snapped point (from another edge)
label pointi = vertsToPoint.lookup(vertices, -1);
if (pointi == -1)
{
- pointi = pointToVerts.size();
+ bool addNewPoint(true);
+ bool edgeIsDiag(pointEdgeType & POINT_DIAGONAL);
+
+ const label snapPointi =
+ (
+ (pointEdgeType & POINT_SNAP0) ? vertices.first()
+ : (pointEdgeType & POINT_SNAP1) ? vertices.second()
+ : -1
+ );
+
+ if (snapPointi == -1)
+ {
+ // No snapped point
+ pointi = pointToVerts.size();
+ pointToVerts.append(vertices);
+ }
+ else
+ {
+ // Snapped point. No corresponding face or diagonal
+ facei = -1;
+ edgeIsDiag = false;
+
+ pointi = snapVertsLookup.lookup(snapPointi, -1);
+ addNewPoint = (pointi == -1);
+ if (addNewPoint)
+ {
+ pointi = pointToVerts.size();
+ snapVertsLookup.insert(snapPointi, pointi);
+ pointToVerts.append(edge(snapPointi, snapPointi));
+ }
+ }
+
+ if (addNewPoint)
+ {
+ pointToFace.append(facei);
+ pointFromDiag.append(edgeIsDiag);
+ }
- pointToVerts.append(vertices);
- pointToFace.append(facei);
- pointFromDiag.append(edgeIsDiag);
vertsToPoint.insert(vertices, pointi);
}
@@ -268,379 +414,370 @@ Foam::label Foam::isoSurfaceTopo::generatePoint
}
-void Foam::isoSurfaceTopo::generateTriPoints
+void Foam::isoSurfaceTopo::generateTetCutPoints
(
const label facei,
const int tetCutIndex,
const tetCell& tetLabels,
- const FixedList<bool, 6>& edgeIsDiag,// Per tet edge whether is face diag
+
+ // Per tet edge whether is face diag etc
+ const FixedList<uint8_t, 6>& edgeDiagType,
DynamicList<edge>& pointToVerts,
DynamicList<label>& pointToFace,
DynamicList<bool>& pointFromDiag,
EdgeMap<label>& vertsToPoint,
+ Map<label>& snapVertsLookup,
+
DynamicList<label>& verts // Every three verts is new triangle
-) const
+)
{
+ bool flip(false);
+
// Form the vertices of the triangles for each case
- switch (tetCutIndex)
+ switch (tetCutIndex & 0x0F)
{
case 0x00:
case 0x0F:
break;
- case 0x01:
- case 0x0E:
+ // Cut point 0
+ case 0x0E: flip = true; [[fallthrough]]; // Point 0 above cut
+ case 0x01: // Point 0 below cut
{
- verts.append
+ const label cutA
(
generatePoint
(
facei,
- edgeIsDiag[0],
- tetLabels.edge(0),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[0] | SNAP_INDEX_DECODE(0, tetCutIndex)),
+ tetLabels.edge(0), // 0 -> 1
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append
+ const label cutB
(
generatePoint
(
facei,
- edgeIsDiag[1],
- tetLabels.edge(1),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[1] | SNAP_INDEX_DECODE(1, tetCutIndex)),
+ tetLabels.edge(1), // 0 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append
+ const label cutC
(
generatePoint
(
facei,
- edgeIsDiag[2],
- tetLabels.edge(2),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[2] | SNAP_INDEX_DECODE(2, tetCutIndex)),
+ tetLabels.edge(2), // 0 -> 3
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- if (tetCutIndex == 0x0E)
- {
- // Flip normals
- const label sz = verts.size();
- std::swap(verts[sz-2], verts[sz-1]);
- }
+ appendTriLabels(verts, cutA, cutB, cutC, flip);
}
break;
- case 0x02:
- case 0x0D:
+ // Cut point 1
+ case 0x0D: flip = true; [[fallthrough]]; // Point 1 above cut
+ case 0x02: // Point 1 below cut
{
- verts.append
+ const label cutA
(
generatePoint
(
facei,
- edgeIsDiag[0],
- tetLabels.reverseEdge(0),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[0] | SNAP_INDEX_DECODE(0, tetCutIndex)),
+ tetLabels.edge(0), // 0 -> 1
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append
+ const label cutB
(
generatePoint
(
facei,
- edgeIsDiag[3],
- tetLabels.reverseEdge(3),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[3] | SNAP_INDEX_DECODE(3, tetCutIndex)),
+ tetLabels.edge(3), // 3 -> 1
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append
+ const label cutC
(
generatePoint
(
facei,
- edgeIsDiag[4],
- tetLabels.edge(4),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[4] | SNAP_INDEX_DECODE(4, tetCutIndex)),
+ tetLabels.edge(4), // 1 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- if (tetCutIndex == 0x0D)
- {
- // Flip normals
- const label sz = verts.size();
- std::swap(verts[sz-2], verts[sz-1]);
- }
+ appendTriLabels(verts, cutA, cutB, cutC, flip);
}
break;
- case 0x03:
- case 0x0C:
+ // Cut point 0/1 | 2/3
+ case 0x0C: flip = true; [[fallthrough]]; // Point 0/1 above cut
+ case 0x03: // Point 0/1 below cut
{
- const label p0p2
+ const label cutA
(
generatePoint
(
facei,
- edgeIsDiag[1],
- tetLabels.edge(1),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[1] | SNAP_INDEX_DECODE(1, tetCutIndex)),
+ tetLabels.edge(1), // 0 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- const label p1p3
+ const label cutB
(
generatePoint
(
facei,
- edgeIsDiag[3],
- tetLabels.reverseEdge(3),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[2] | SNAP_INDEX_DECODE(2, tetCutIndex)),
+ tetLabels.edge(2), // 0 -> 3
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
-
- verts.append
+ const label cutC
(
generatePoint
(
facei,
- edgeIsDiag[2],
- tetLabels.edge(2),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[3] | SNAP_INDEX_DECODE(3, tetCutIndex)),
+ tetLabels.edge(3), // 3 -> 1
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append(p1p3);
- verts.append(p0p2);
- verts.append(p1p3);
- verts.append
+ const label cutD
(
generatePoint
(
facei,
- edgeIsDiag[4],
- tetLabels.edge(4),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[4] | SNAP_INDEX_DECODE(4, tetCutIndex)),
+ tetLabels.edge(4), // 1 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append(p0p2);
- if (tetCutIndex == 0x0C)
- {
- // Flip normals
- const label sz = verts.size();
- std::swap(verts[sz-5], verts[sz-4]);
- std::swap(verts[sz-2], verts[sz-1]);
- }
+ appendTriLabels(verts, cutA, cutB, cutC, flip);
+ appendTriLabels(verts, cutA, cutC, cutD, flip);
}
break;
- case 0x04:
- case 0x0B:
+ // Cut point 2
+ case 0x0B: flip = true; [[fallthrough]]; // Point 2 above cut
+ case 0x04: // Point 2 below cut
{
- verts.append
+ const label cutA
(
generatePoint
(
facei,
- edgeIsDiag[1],
- tetLabels.reverseEdge(1),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[1] | SNAP_INDEX_DECODE(1, tetCutIndex)),
+ tetLabels.edge(1), // 0 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append
+ const label cutB
(
generatePoint
(
facei,
- edgeIsDiag[4],
- tetLabels.reverseEdge(4),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[4] | SNAP_INDEX_DECODE(4, tetCutIndex)),
+ tetLabels.edge(4), // 1 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append
+ const label cutC
(
generatePoint
(
facei,
- edgeIsDiag[5],
- tetLabels.reverseEdge(5),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[5] | SNAP_INDEX_DECODE(5, tetCutIndex)),
+ tetLabels.edge(5), // 3 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- if (tetCutIndex == 0x0B)
- {
- // Flip normals
- const label sz = verts.size();
- std::swap(verts[sz-2], verts[sz-1]);
- }
+ appendTriLabels(verts, cutA, cutB, cutC, flip);
}
break;
- case 0x05:
- case 0x0A:
+ // Cut point 0/2 | 1/3
+ case 0x0A: flip = true; [[fallthrough]]; // Point 0/2 above cut
+ case 0x05: // Point 0/2 below cut
{
- const label p0p1
+ const label cutA
(
generatePoint
(
facei,
- edgeIsDiag[0],
- tetLabels.edge(0),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[0] | SNAP_INDEX_DECODE(0, tetCutIndex)),
+ tetLabels.edge(0), // 0 -> 1
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- const label p2p3
+ const label cutB
(
generatePoint
(
facei,
- edgeIsDiag[5],
- tetLabels.reverseEdge(5),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[4] | SNAP_INDEX_DECODE(4, tetCutIndex)),
+ tetLabels.edge(4), // 1 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
-
- verts.append(p0p1);
- verts.append(p2p3);
- verts.append
+ const label cutC
(
generatePoint
(
facei,
- edgeIsDiag[2],
- tetLabels.edge(2),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[5] | SNAP_INDEX_DECODE(5, tetCutIndex)),
+ tetLabels.edge(5), // 3 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append(p0p1);
- verts.append
+ const label cutD
(
generatePoint
(
facei,
- edgeIsDiag[4],
- tetLabels.edge(4),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[2] | SNAP_INDEX_DECODE(2, tetCutIndex)),
+ tetLabels.edge(2), // 0 -> 3
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append(p2p3);
- if (tetCutIndex == 0x0A)
- {
- // Flip normals
- const label sz = verts.size();
- std::swap(verts[sz-5], verts[sz-4]);
- std::swap(verts[sz-2], verts[sz-1]);
- }
+ appendTriLabels(verts, cutA, cutB, cutC, flip);
+ appendTriLabels(verts, cutA, cutC, cutD, flip);
}
break;
- case 0x06:
- case 0x09:
+ // Cut point 1/2 | 0/3
+ case 0x09: flip = true; [[fallthrough]]; // Point 1/2 above cut
+ case 0x06: // Point 1/2 below cut
{
- const label p0p1
+ const label cutA
(
generatePoint
(
facei,
- edgeIsDiag[0],
- tetLabels.edge(0),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[0] | SNAP_INDEX_DECODE(0, tetCutIndex)),
+ tetLabels.edge(0), // 0 -> 1
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- const label p2p3
+ const label cutB
(
generatePoint
(
facei,
- edgeIsDiag[5],
- tetLabels.reverseEdge(5),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[3] | SNAP_INDEX_DECODE(3, tetCutIndex)),
+ tetLabels.edge(3), // 3 -> 1
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
-
- verts.append(p0p1);
- verts.append
+ const label cutC
(
generatePoint
(
facei,
- edgeIsDiag[3],
- tetLabels.reverseEdge(3),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[5] | SNAP_INDEX_DECODE(5, tetCutIndex)),
+ tetLabels.edge(5), // 3 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append(p2p3);
- verts.append(p0p1);
- verts.append(p2p3);
- verts.append
+ const label cutD
(
generatePoint
(
facei,
- edgeIsDiag[1],
- tetLabels.edge(1),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[1] | SNAP_INDEX_DECODE(1, tetCutIndex)),
+ tetLabels.edge(1), // 0 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- if (tetCutIndex == 0x09)
- {
- // Flip normals
- const label sz = verts.size();
- std::swap(verts[sz-5], verts[sz-4]);
- std::swap(verts[sz-2], verts[sz-1]);
- }
+ appendTriLabels(verts, cutA, cutB, cutC, flip);
+ appendTriLabels(verts, cutA, cutC, cutD, flip);
}
break;
- case 0x08:
- case 0x07:
+ // Cut point 3
+ case 0x07: flip = true; [[fallthrough]]; // Point 3 above cut
+ case 0x08: // Point 3 below cut
{
- verts.append
+ const label cutA
(
generatePoint
(
facei,
- edgeIsDiag[2],
- tetLabels.reverseEdge(2),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[2] | SNAP_INDEX_DECODE(2, tetCutIndex)),
+ tetLabels.edge(2), // 0 -> 3
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append
+ const label cutB
(
generatePoint
(
facei,
- edgeIsDiag[5],
- tetLabels.reverseEdge(5),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[5] | SNAP_INDEX_DECODE(5, tetCutIndex)),
+ tetLabels.edge(5), // 3 -> 2
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- verts.append
+ const label cutC
(
generatePoint
(
facei,
- edgeIsDiag[3],
- tetLabels.edge(3),
- pointToVerts, pointToFace, pointFromDiag, vertsToPoint
+ (edgeDiagType[3] | SNAP_INDEX_DECODE(3, tetCutIndex)),
+ tetLabels.edge(3), // 3 -> 1
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup
)
);
- if (tetCutIndex == 0x07)
- {
- // Flip normals
- const label sz = verts.size();
- std::swap(verts[sz-2], verts[sz-1]);
- }
+
+ appendTriLabels(verts, cutA, cutB, cutC, flip);
}
break;
}
}
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+// Requires mesh_, tetBasePtIs_
void Foam::isoSurfaceTopo::generateTriPoints
(
const label celli,
@@ -651,21 +788,21 @@ void Foam::isoSurfaceTopo::generateTriPoints
DynamicList<bool>& pointFromDiag,
EdgeMap<label>& vertsToPoint,
- DynamicList<label>& verts,
- DynamicList<label>& faceLabels
+ Map<label>& snapVertsLookup,
+
+ DynamicList<label>& verts
) const
{
const faceList& faces = mesh_.faces();
const labelList& faceOwner = mesh_.faceOwner();
const cell& cFaces = mesh_.cells()[celli];
+ const bool doSnap = this->snap();
if (isTet)
{
// For tets don't do cell-centre decomposition, just use the
// tet points and values
- const label startTrii = verts.size();
-
const label facei = cFaces[0];
const face& f0 = faces[facei];
@@ -682,7 +819,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
}
}
- label p0 = f0[0];
+ const label p0 = f0[0];
label p1 = f0[1];
label p2 = f0[2];
@@ -691,7 +828,7 @@ void Foam::isoSurfaceTopo::generateTriPoints
std::swap(p1, p2);
}
- generateTriPoints
+ generateTetCutPoints
(
facei,
getTetCutIndex
@@ -700,22 +837,17 @@ void Foam::isoSurfaceTopo::generateTriPoints
pVals_[p1],
pVals_[p2],
pVals_[oppositeI],
- iso_
+ iso_,
+ doSnap
),
tetCell(p0, p1, p2, oppositeI),
- FixedList<bool, 6>(false), // edgeIsDiag = false
+ FixedList<uint8_t, 6>(POINT_NONE), // edgeDiagType
+
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup,
- pointToVerts,
- pointToFace,
- pointFromDiag,
- vertsToPoint,
verts // Every three verts is new triangle
);
-
- for (label nTris = (verts.size()-startTrii)/3; nTris; --nTris)
- {
- faceLabels.append(facei);
- }
}
else
{
@@ -730,8 +862,6 @@ void Foam::isoSurfaceTopo::generateTriPoints
continue;
}
- const label startTrii = verts.size();
-
const face& f = faces[facei];
label fp0 = tetBasePtIs_[facei];
@@ -742,29 +872,29 @@ void Foam::isoSurfaceTopo::generateTriPoints
fp0 = 0;
}
+ const label p0 = f[fp0];
label fp = f.fcIndex(fp0);
for (label i = 2; i < f.size(); ++i)
{
const label nextFp = f.fcIndex(fp);
- FixedList<bool, 6> edgeIsDiag(false);
+ FixedList<uint8_t, 6> edgeDiagType(POINT_NONE);
- label p0 = f[fp0];
label p1 = f[fp];
label p2 = f[nextFp];
if (faceOwner[facei] == celli)
{
std::swap(p1, p2);
- if (i != 2) edgeIsDiag[1] = true;
- if (i != f.size()-1) edgeIsDiag[0] = true;
+ if (i != 2) edgeDiagType[1] = POINT_DIAGONAL;
+ if (i != f.size()-1) edgeDiagType[0] = POINT_DIAGONAL;
}
else
{
- if (i != 2) edgeIsDiag[0] = true;
- if (i != f.size()-1) edgeIsDiag[1] = true;
+ if (i != 2) edgeDiagType[0] = POINT_DIAGONAL;
+ if (i != f.size()-1) edgeDiagType[1] = POINT_DIAGONAL;
}
- generateTriPoints
+ generateTetCutPoints
(
facei,
getTetCutIndex
@@ -773,25 +903,20 @@ void Foam::isoSurfaceTopo::generateTriPoints
pVals_[p1],
pVals_[p2],
cVals_[celli],
- iso_
+ iso_,
+ doSnap
),
tetCell(p0, p1, p2, mesh_.nPoints()+celli),
- edgeIsDiag,
+ edgeDiagType,
+
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup,
- pointToVerts,
- pointToFace,
- pointFromDiag,
- vertsToPoint,
verts // Every three verts is new triangle
);
fp = nextFp;
}
-
- for (label nTris = (verts.size()-startTrii)/3; nTris; --nTris)
- {
- faceLabels.append(facei);
- }
}
}
}
@@ -807,7 +932,7 @@ void Foam::isoSurfaceTopo::triangulateOutside
const labelUList& pointToFace,
const label cellID,
- // outputs
+ // output
DynamicList<face>& compactFaces,
DynamicList<label>& compactCellIDs
)
@@ -916,7 +1041,6 @@ void Foam::isoSurfaceTopo::removeInsidePoints
pp,
pointFromDiag,
pointToFace,
- //protectedFace,
celli,
compactFaces,
@@ -1035,7 +1159,7 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
false
),
fvmesh,
- dimensionedScalar(dimless, Zero)
+ dimensionedScalar(dimless)
);
auto& debugFld = debugField.primitiveFieldRef();
@@ -1045,36 +1169,39 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
debugFld[celli] = cellCutType_[celli];
}
- Pout<< "Writing cut types:"
- << debugField.objectPath() << endl;
-
+ Info<< "Writing cut types: " << debugField.objectRelPath() << nl;
debugField.write();
}
// Per cell: 5 pyramids cut, each generating 2 triangles
- // - pointToVerts : from generated iso point to originating mesh verts
+
+ // From generated iso point to originating mesh verts
DynamicList<edge> pointToVerts(10*nCutCells);
- // - pointToFace : from generated iso point to originating mesh face
+
+ // From generated iso point to originating mesh face, -1 for snapped
DynamicList<label> pointToFace(10*nCutCells);
- // - pointFromDiag : if generated iso point is on face diagonal
+
+ // Generated iso point is on face diagonal
DynamicList<bool> pointFromDiag(10*nCutCells);
// Per cell: number of intersected edges:
// - four faces cut so 4 mesh edges + 4 face-diagonal edges
// - 4 of the pyramid edges
EdgeMap<label> vertsToPoint(12*nCutCells);
- DynamicList<label> verts(12*nCutCells);
- // Per cell: 5 pyramids cut (since only one pyramid not cut)
- DynamicList<label> faceLabels(5*nCutCells);
- DynamicList<label> cellLabels(5*nCutCells);
+ Map<label> snapVertsLookup(0);
+ if (this->snap())
+ {
+ // Some, but not all, cells may have point snapping
+ snapVertsLookup.resize(4*nCutCells);
+ }
+ DynamicList<label> verts(12*nCutCells);
labelList startTri(mesh_.nCells()+1, Zero);
-
for (label celli = 0; celli < mesh_.nCells(); ++celli)
{
- startTri[celli] = faceLabels.size();
+ startTri[celli] = (verts.size()/3);
if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
{
generateTriPoints
@@ -1084,27 +1211,19 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
// Same as tetMatcher::test(mesh_, celli),
bool(cellCutType_[celli] & cutType::TETCUT), // isTet
- pointToVerts,
- pointToFace,
- pointFromDiag,
+ pointToVerts, pointToFace, pointFromDiag,
+ vertsToPoint, snapVertsLookup,
- vertsToPoint,
- verts,
- faceLabels
+ verts
);
-
- for (label i = startTri[celli]; i < faceLabels.size(); ++i)
- {
- cellLabels.append(celli);
- }
}
}
- startTri[mesh_.nCells()] = faceLabels.size();
+ startTri.last() = (verts.size()/3);
+ // Not required after this stage
+ tetBasePtIs_.clear();
pointToVerts_.transfer(pointToVerts);
- meshCells_.transfer(cellLabels);
- pointToFace_.transfer(pointToFace);
pointField allPoints
(
@@ -1116,26 +1235,42 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
);
- // Assign to MeshedSurface
- faceList allTris(faceLabels.size());
- label verti = 0;
- for (face& allTri : allTris)
+ // List of faces, and cells cut by them
+ meshCells_.resize(startTri.last());
+ for (label celli = 0; celli < startTri.size()-1; ++celli)
{
- allTri.resize(3);
- allTri[0] = verts[verti++];
- allTri[1] = verts[verti++];
- allTri[2] = verts[verti++];
+ // All triangles for the current cell
+ SubList<label>
+ (
+ meshCells_,
+ (startTri[celli+1] - startTri[celli]),
+ startTri[celli]
+ ) = celli;
}
- surfZoneList allZones(one{}, surfZone("allFaces", allTris.size()));
+ // From list of vertices to triangulate faces
+ faceList allTris(startTri.last());
+ {
+ label verti = 0;
+ for (face& f : allTris)
+ {
+ f.resize(3);
+ f[0] = verts[verti++];
+ f[1] = verts[verti++];
+ f[2] = verts[verti++];
+ }
+ verts.clearStorage();
+ }
+
+ // Assign to MeshedSurface
Mesh::clear();
Mesh updated
(
std::move(allPoints),
std::move(allTris),
- std::move(allZones)
+ surfZoneList(one{}, surfZone("allFaces", allTris.size()))
);
Mesh::transfer(updated);
@@ -1143,7 +1278,7 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
// - generated faces and points are assigned to *this
// - per point we know:
// - pointOnDiag: whether it is on a face-diagonal edge
- // - pointToFace_: from what pyramid (cell+face) it was produced
+ // - pointToFace: from what pyramid (cell+face) it was produced
// (note that the pyramid faces are shared between multiple mesh faces)
// - pointToVerts_ : originating mesh vertex or cell centre
@@ -1169,14 +1304,13 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
(params.filter() == filterType::DIAGCELL ? true : false),
pointFromDiag,
- pointToFace_,
+ pointToFace,
startTri,
pointCompactMap,
compactCellIDs
);
pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointCompactMap)();
- pointToFace_ = UIndirectList<label>(pointToFace_, pointCompactMap)();
meshCells_.transfer(compactCellIDs);
pointCompactMap.clearStorage();
@@ -1194,6 +1328,7 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
// Not required after this stage
pointFromDiag.clearStorage();
+ pointToFace.clearStorage();
if (params.filter() == filterType::DIAGCELL)
@@ -1270,10 +1405,9 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
const labelList& eFaces = edgeFaces[edgei];
if (eFaces.size() == 1)
{
- // Open edge. Check that vertices do not originate
- // from a boundary face
- const edge& e = surfEdges[edgei];
+ // Open edge (not originating from a boundary face)
+ const edge& e = surfEdges[edgei];
const edge& verts0 = pointToVerts_[mp[e.first()]];
const edge& verts1 = pointToVerts_[mp[e.second()]];
@@ -1344,7 +1478,6 @@ void Foam::isoSurfaceTopo::inplaceSubsetMesh(const bitSet& include)
meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
- pointToFace_ = UIndirectList<label>(pointToFace_, pointMap)();
}
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.H b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
index 1a5bba527cd..5c5d72ff360 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2019 OpenFOAM Foundation
- Copyright (C) 2019-2020 OpenCFD Ltd.
+ Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -60,17 +60,29 @@ class isoSurfaceTopo
:
public isoSurfaceBase
{
+ // Data Types
+
+ //- Masks for point type on edge
+ enum pointEdgeMask : uint8_t
+ {
+ POINT_NONE = 0, //!< No information
+ POINT_SNAP0 = 1, //!< Point is snap to first
+ POINT_SNAP1 = 2, //!< Point is snap to second
+ POINT_DIAGONAL = 4 //!< Point is on face diagonal
+ };
+
+
// Private Data
- //- Corrected version of tetBasePtIs
+ //- Corrected version of tetBasePtIs.
+ // Used during construction
labelList tetBasePtIs_;
- //- Per point: originating mesh vertex/cc. See encoding above
+ //- Per point: originating mesh vertex/cell centre combination.
+ // Vertices less than nPoints are mesh vertices,
+ // duplicate vertices indicate a snapped point
edgeList pointToVerts_;
- //- For every point the originating face in mesh
- labelList pointToFace_;
-
//- The cell cut type
List<cutType> cellCutType_;
@@ -86,47 +98,60 @@ class isoSurfaceTopo
void fixTetBasePtIs();
//- Generate single point on edge
- label generatePoint
+ static label generatePoint
(
- const label facei,
- const bool edgeIsDiag,
+ label facei, //!< Originating mesh face for the iso-point
+ const int pointEdgeType, // Snap to 0/1, diagonal...
const edge& vertices,
+ // Output
DynamicList<edge>& pointToVerts,
DynamicList<label>& pointToFace,
DynamicList<bool>& pointFromDiag,
- EdgeMap<label>& vertsToPoint
- ) const;
- //- Generate triangles from tet
- void generateTriPoints
+ // Bookkeeping
+ EdgeMap<label>& vertsToPoint,
+ Map<label>& snapVertsLookup
+ );
+
+ //- Generate triangles from tet cut
+ static void generateTetCutPoints
(
const label facei,
- const int tetCutIndex, //!< Encoded tet cuts. getTetCutIndex()
+ const int tetCutIndex, //!< Encoded tet cuts
const tetCell& tetLabels,
- const FixedList<bool, 6>& edgeIsDiag,
+ const FixedList<uint8_t, 6>& edgeDiagType,
+ // Output
DynamicList<edge>& pointToVerts,
DynamicList<label>& pointToFace,
DynamicList<bool>& pointFromDiag,
+ // Bookkeeping
EdgeMap<label>& vertsToPoint,
+ Map<label>& snapVertsLookup,
+
+ // Output
DynamicList<label>& verts
- ) const;
+ );
- //- Generate triangles from cell
+ //- Generate triangles from cell.
void generateTriPoints
(
const label celli,
const bool isTet,
+ // Output
DynamicList<edge>& pointToVerts,
DynamicList<label>& pointToFace,
DynamicList<bool>& pointFromDiag,
+ // Bookkeeping
EdgeMap<label>& vertsToPoint,
- DynamicList<label>& verts,
- DynamicList<label>& faceLabels
+ Map<label>& snapVertsLookup,
+
+ // Output
+ DynamicList<label>& verts
) const;
@@ -173,12 +198,12 @@ protected:
// Sampling
- //- Interpolates cCoords,pCoords.
+ //- Interpolates cellData and pointData fields
template<class Type>
tmp<Field<Type>> interpolateTemplate
(
- const Field<Type>& cCoords,
- const Field<Type>& pCoords
+ const Field<Type>& cellData,
+ const Field<Type>& pointData
) const;
public:
@@ -213,14 +238,10 @@ public:
// Member Functions
- //- For every point originating face (pyramid) in mesh
- const labelList& pointToFace() const
- {
- return pointToFace_;
- }
-
- //- Per point: originating mesh vertex/cc. See encoding above
- const edgeList& pointToVerts() const
+ //- Per point: originating mesh vertex/cell centre combination.
+ // Vertices less than nPoints are mesh vertices (encoding above),
+ // duplicate vertices indicate a snapped point
+ const edgeList& pointToVerts() const noexcept
{
return pointToVerts_;
}
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
index 8ffec74e8db..5056fb38a5e 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopoTemplates.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2018 OpenFOAM Foundation
- Copyright (C) 2020 OpenCFD Ltd.
+ Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -32,8 +32,8 @@ template<class Type>
Foam::tmp<Foam::Field<Type>>
Foam::isoSurfaceTopo::interpolateTemplate
(
- const Field<Type>& cellCoords,
- const Field<Type>& pointCoords
+ const Field<Type>& cellData,
+ const Field<Type>& pointData
) const
{
auto tfld = tmp<Field<Type>>::New(pointToVerts_.size());
@@ -41,41 +41,50 @@ Foam::isoSurfaceTopo::interpolateTemplate
forAll(pointToVerts_, i)
{
+ const edge& verts = pointToVerts_[i];
+ Type& val = fld[i];
+
scalar s0;
- Type p0;
+ Type v0;
{
- label idx = pointToVerts_[i].first();
+ label idx = verts.first();
if (idx < mesh_.nPoints())
{
// Point index
s0 = pVals_[idx];
- p0 = pointCoords[idx];
+ v0 = pointData[idx];
}
else
{
// Cell index
idx -= mesh_.nPoints();
s0 = cVals_[idx];
- p0 = cellCoords[idx];
+ v0 = cellData[idx];
}
}
scalar s1;
- Type p1;
+ Type v1;
{
- label idx = pointToVerts_[i].second();
- if (idx < mesh_.nPoints())
+ label idx = verts.second();
+ if (idx == verts.first())
+ {
+ // Duplicate index (ie, snapped)
+ val = v0;
+ continue;
+ }
+ else if (idx < mesh_.nPoints())
{
// Point index
s1 = pVals_[idx];
- p1 = pointCoords[idx];
+ v1 = pointData[idx];
}
else
{
// Cell index
idx -= mesh_.nPoints();
s1 = cVals_[idx];
- p1 = cellCoords[idx];
+ v1 = cellData[idx];
}
}
@@ -83,11 +92,11 @@ Foam::isoSurfaceTopo::interpolateTemplate
if (mag(d) > VSMALL)
{
const scalar s = (iso_-s0)/d;
- fld[i] = s*p1+(1.0-s)*p0;
+ val = s*v1+(1.0-s)*v0;
}
else
{
- fld[i] = 0.5*(p0+p1);
+ val = 0.5*(v0+v1);
}
}
--
GitLab