diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.C b/src/sampling/surface/isoSurface/isoSurfaceCell.C
index 4e27f0084e84e77770abda930d939af1bcc4c0c2..9ac005c1bbb68775928cff43c5fac4f51f34c5b3 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceCell.C
@@ -45,6 +45,29 @@ namespace Foam
}
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ // Does any edge of triangle cross iso value?
+ inline static bool isTriCut
+ (
+ const label a,
+ const label b,
+ const label c,
+ const scalarField& pointValues,
+ const scalar isoValue
+ )
+ {
+ const bool aLower = (pointValues[a] < isoValue);
+ const bool bLower = (pointValues[b] < isoValue);
+ const bool cLower = (pointValues[c] < isoValue);
+
+ return !(aLower == bLower && aLower == cLower);
+ }
+}
+
+
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::scalar Foam::isoSurfaceCell::isoFraction
@@ -64,20 +87,6 @@ Foam::scalar Foam::isoSurfaceCell::isoFraction
}
-bool Foam::isoSurfaceCell::isTriCut
-(
- const triFace& tri,
- const scalarField& pointValues
-) const
-{
- const bool aLower = (pointValues[tri[0]] < iso_);
- const bool bLower = (pointValues[tri[1]] < iso_);
- const bool cLower = (pointValues[tri[2]] < iso_);
-
- return !(aLower == bLower && aLower == cLower);
-}
-
-
Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
(
const bitSet& isTet,
@@ -101,9 +110,7 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
for (label fp = 1; fp < f.size() - 1; ++fp)
{
- triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-
- if (isTriCut(tri, pointValues))
+ if (isTriCut(f[0], f[fp], f[f.fcIndex(fp)], pointValues, iso_))
{
return CUT;
}
@@ -112,6 +119,7 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
return NOTCUT;
}
+
const bool cellLower = (cellValues[celli] < iso_);
// First check if there is any cut in cell
@@ -137,13 +145,16 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
}
// Check (triangulated) face edges
- const label fp0 = mesh_.tetBasePtIs()[facei];
+ // With fallback for problem decompositions
+ const label fp0 =
+ (mesh_.tetBasePtIs()[facei] < 0 ? 0 : mesh_.tetBasePtIs()[facei]);
+
label fp = f.fcIndex(fp0);
for (label i = 2; i < f.size(); ++i)
{
const label nextFp = f.fcIndex(fp);
- if (isTriCut(triFace(f[fp0], f[fp], f[nextFp]), pointValues))
+ if (isTriCut(f[fp0], f[fp], f[nextFp], pointValues, iso_))
{
edgeCut = true;
break;
@@ -167,21 +178,21 @@ Foam::isoSurfaceCell::cellCutType Foam::isoSurfaceCell::calcCutType
const labelList& cPoints = mesh_.cellPoints(celli);
- label nPyrEdgeCuts = 0;
+ label nPyrCuts = 0;
for (const label pointi : cPoints)
{
- if (cellLower != (pointValues[pointi] < iso_))
+ if ((pointValues[pointi] < iso_) != cellLower)
{
- ++nPyrEdgeCuts;
+ ++nPyrCuts;
}
}
- if (nPyrEdgeCuts == cPoints.size())
+ if (nPyrCuts == cPoints.size())
{
return SPHERE;
}
- else if (nPyrEdgeCuts)
+ else if (nPyrCuts)
{
// There is a pyramid edge cut. E.g. lopping off a tet from a corner
return CUT;
@@ -206,7 +217,7 @@ void Foam::isoSurfaceCell::calcCutTypes
// Do now to ensure things stay synchronized.
(void)mesh_.tetBasePtIs();
- forAll(mesh_.cells(), celli)
+ forAll(cellCutType_, celli)
{
cellCutType_[celli] = calcCutType(isTet, cVals, pVals, celli);
diff --git a/src/sampling/surface/isoSurface/isoSurfaceCell.H b/src/sampling/surface/isoSurface/isoSurfaceCell.H
index 34e022f5f6588ff12cb00b84f13400e4fb2396f4..b4e4c5076e6e3722002b300ea44c3c60e18638bb 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceCell.H
@@ -128,12 +128,6 @@ class isoSurfaceCell
const scalar s1
) const;
- //- Does any edge of triangle cross iso value?
- bool isTriCut
- (
- const triFace& tri,
- const scalarField& pointValues
- ) const;
//- Determine whether cell is cut
cellCutType calcCutType
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.C b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
index 5d52e79150b26b12df936ddfd235ce46a0e293a7..744d995c7b0e9c9343b6b12f47daa4b98b16811b 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.C
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.C
@@ -42,22 +42,31 @@ namespace Foam
}
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
-bool Foam::isoSurfaceTopo::isTriCut
-(
- const triFace& tri,
- const scalarField& pointValues
-) const
+namespace Foam
{
- const bool aLower = (pointValues[tri[0]] < iso_);
- const bool bLower = (pointValues[tri[1]] < iso_);
- const bool cLower = (pointValues[tri[2]] < iso_);
+ // Does any edge of triangle cross iso value?
+ inline static bool isTriCut
+ (
+ const label a,
+ const label b,
+ const label c,
+ const scalarField& pointValues,
+ const scalar isoValue
+ )
+ {
+ const bool aLower = (pointValues[a] < isoValue);
+ const bool bLower = (pointValues[b] < isoValue);
+ const bool cLower = (pointValues[c] < isoValue);
- return !(aLower == bLower && aLower == cLower);
+ return !(aLower == bLower && aLower == cLower);
+ }
}
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
Foam::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
(
const bool isTet,
@@ -79,9 +88,7 @@ Foam::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
for (label fp = 1; fp < f.size() - 1; ++fp)
{
- const triFace tri(f[0], f[fp], f[f.fcIndex(fp)]);
-
- if (isTriCut(tri, pVals_))
+ if (isTriCut(f[0], f[fp], f[f.fcIndex(fp)], pVals_, iso_))
{
return CUT;
}
@@ -89,89 +96,85 @@ Foam::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
}
return NOTCUT;
}
- else
- {
- const bool cellLower = (cVals_[celli] < iso_);
- // First check if there is any cut in cell
- bool edgeCut = false;
- for (const label facei : cFaces)
- {
- const face& f = mesh_.faces()[facei];
-
- // Check pyramids cut
- for (const label pointi : f)
- {
- if ((pVals_[pointi] < iso_) != cellLower)
- {
- edgeCut = true;
- break;
- }
- }
+ const bool cellLower = (cVals_[celli] < iso_);
- if (edgeCut)
- {
- break;
- }
+ // First check if there is any cut in cell
+ bool edgeCut = false;
- label fp0 = tetBasePtIs_[facei];
+ for (const label facei : cFaces)
+ {
+ const face& f = mesh_.faces()[facei];
- // Fall back for problem decompositions
- if (fp0 < 0)
+ // Check pyramids cut
+ for (const label pointi : f)
+ {
+ if ((pVals_[pointi] < iso_) != cellLower)
{
- fp0 = 0;
+ edgeCut = true;
+ break;
}
+ }
- label fp = f.fcIndex(fp0);
- for (label i = 2; i < f.size(); ++i)
- {
- label nextFp = f.fcIndex(fp);
+ if (edgeCut)
+ {
+ break;
+ }
- if (isTriCut(triFace(f[fp0], f[fp], f[nextFp]), pVals_))
- {
- edgeCut = true;
- break;
- }
+ // Check (triangulated) face edges
+ // With fallback for problem decompositions
+ const label fp0 = (tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei]);
- fp = nextFp;
- }
+ label fp = f.fcIndex(fp0);
+ for (label i = 2; i < f.size(); ++i)
+ {
+ const label nextFp = f.fcIndex(fp);
- if (edgeCut)
+ if (isTriCut(f[fp0], f[fp], f[nextFp], pVals_, iso_))
{
+ edgeCut = true;
break;
}
+
+ fp = nextFp;
}
if (edgeCut)
{
- // Count actual cuts (expensive since addressing needed)
- // Note: not needed if you don't want to preserve maxima/minima
- // centred around cellcentre. In that case just always return CUT
+ break;
+ }
+ }
- const labelList& cPoints = mesh_.cellPoints(celli);
+ if (edgeCut)
+ {
+ // Count actual cuts (expensive since addressing needed)
+ // Note: not needed if you don't want to preserve maxima/minima
+ // centred around cellcentre. In that case just always return CUT
- label nPyrEdgeCuts = 0;
+ const labelList& cPoints = mesh_.cellPoints(celli);
- for (const label pointi : cPoints)
- {
- if ((pVals_[pointi] < iso_) != cellLower)
- {
- ++nPyrEdgeCuts;
- }
- }
+ label nPyrCuts = 0;
- if (nPyrEdgeCuts == cPoints.size())
- {
- return SPHERE;
- }
- else if (nPyrEdgeCuts)
+ for (const label pointi : cPoints)
+ {
+ if ((pVals_[pointi] < iso_) != cellLower)
{
- return CUT;
+ ++nPyrCuts;
}
}
- return NOTCUT;
+
+ if (nPyrCuts == cPoints.size())
+ {
+ return SPHERE;
+ }
+ else if (nPyrCuts)
+ {
+ return CUT;
+ }
}
+
+ return NOTCUT;
}
@@ -181,11 +184,10 @@ Foam::label Foam::isoSurfaceTopo::calcCutTypes
List<cellCutType>& cellCutTypes
)
{
- const cellList& cells = mesh_.cells();
+ cellCutTypes.setSize(mesh_.nCells());
- cellCutTypes.setSize(cells.size());
label nCutCells = 0;
- forAll(cells, celli)
+ forAll(cellCutTypes, celli)
{
cellCutTypes[celli] = calcCutType(tet.isA(mesh_, celli), celli);
@@ -394,22 +396,20 @@ Foam::label Foam::isoSurfaceTopo::generatePoint
EdgeMap<label>& vertsToPoint
) const
{
- EdgeMap<label>::const_iterator edgeFnd = vertsToPoint.find(vertices);
- if (edgeFnd != vertsToPoint.end())
- {
- return edgeFnd();
- }
- else
+ const auto edgeFnd = vertsToPoint.cfind(vertices);
+ if (edgeFnd.found())
{
- // Generate new point
- label pointi = pointToVerts.size();
-
- pointToVerts.append(vertices);
- pointToFace.append(facei);
- pointFromDiag.append(edgeIsDiag);
- vertsToPoint.insert(vertices, pointi);
- return pointi;
+ return edgeFnd.val();
}
+
+ // Generate new point
+ label pointi = pointToVerts.size();
+
+ pointToVerts.append(vertices);
+ pointToFace.append(facei);
+ pointFromDiag.append(edgeIsDiag);
+ vertsToPoint.insert(vertices, pointi);
+ return pointi;
}
@@ -1181,7 +1181,8 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
<< min(pVals_) << " / "
<< max(pVals_) << nl
<< " isoValue : " << iso << nl
- << " filter : " << int(filter) << nl
+ << " filter : " << isoSurfaceTopo::filterNames[filter]
+ << nl
<< " mesh span : " << mesh.bounds().mag() << nl
<< " ignoreCells : " << ignoreCells_.count()
<< " / " << cVals_.size() << nl
@@ -1329,11 +1330,8 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
}
- if (filter == filterType::DIAGCELL && ignoreCells_.empty())
+ if (filter == filterType::DIAGCELL)
{
- // Note that the following only works without cell subsets
- // thus we skip this when ignoreCells_ is non-empty
-
// We remove verts on face diagonals. This is in fact just
// straightening the edges of the face through the cell. This can
// close off 'pockets' of triangles and create open or
@@ -1356,43 +1354,71 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
}
+ // Include faces that would be exposed from mesh subset
+ if (!ignoreCells_.empty())
+ {
+ const labelList& faceOwner = mesh_.faceOwner();
+ const labelList& faceNeighbour = mesh_.faceNeighbour();
+
+ label nMore = 0;
+ for
+ (
+ label facei = 0;
+ facei < mesh.nInternalFaces();
+ ++facei
+ )
+ {
+ int n = 0;
+ if (ignoreCells_.test(faceOwner[facei])) ++n;
+ if (ignoreCells_.test(faceNeighbour[facei])) ++n;
+
+ // Only one of the cells is being ignored.
+ // That means it is an exposed subMesh face.
+ if (n == 1)
+ {
+ isBoundaryPoint.set(mesh.faces()[facei]);
+
+ ++nMore;
+ }
+ }
+ }
+
+
+ bitSet faceSelection;
+
while (true)
{
+ faceSelection.clear();
+ faceSelection.resize(this->size());
+
const labelList& mp = meshPoints();
- bitSet removeFace(this->size());
- label nFaces = 0;
+ const labelListList& edgeFaces = MeshStorage::edgeFaces();
+
+ forAll(edgeFaces, edgei)
{
- const labelListList& edgeFaces =
- MeshStorage::edgeFaces();
- forAll(edgeFaces, edgei)
+ const labelList& eFaces = edgeFaces[edgei];
+ if (eFaces.size() == 1)
{
- const labelList& eFaces = edgeFaces[edgei];
- if (eFaces.size() == 1)
+ // Open edge. Check that vertices do not originate
+ // from a boundary face
+ const edge& e = edges()[edgei];
+ const edge& verts0 = pointToVerts_[mp[e[0]]];
+ const edge& verts1 = pointToVerts_[mp[e[1]]];
+ if
+ (
+ isBoundaryPoint[verts0[0]]
+ && isBoundaryPoint[verts0[1]]
+ && isBoundaryPoint[verts1[0]]
+ && isBoundaryPoint[verts1[1]]
+ )
{
- // Open edge. Check that vertices do not originate
- // from a boundary face
- const edge& e = edges()[edgei];
- const edge& verts0 = pointToVerts_[mp[e[0]]];
- const edge& verts1 = pointToVerts_[mp[e[1]]];
- if
- (
- isBoundaryPoint[verts0[0]]
- && isBoundaryPoint[verts0[1]]
- && isBoundaryPoint[verts1[0]]
- && isBoundaryPoint[verts1[1]]
- )
- {
- // Open edge on boundary face. Keep
- }
- else
- {
- // Open edge. Mark for erosion
- if (removeFace.set(eFaces[0]))
- {
- ++nFaces;
- }
- }
+ // Open edge on boundary face. Keep
+ }
+ else
+ {
+ // Open edge. Mark for erosion
+ faceSelection.set(eFaces[0]);
}
}
}
@@ -1400,24 +1426,18 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
if (debug)
{
Pout<< "isoSurfaceTopo :"
- << " removing " << nFaces
+ << " removing " << faceSelection.count()
<< " faces since on open edges" << endl;
}
- if (returnReduce(nFaces, sumOp<label>()) == 0)
+ if (returnReduce(faceSelection.none(), andOp<bool>()))
{
break;
}
- // Remove the faces
- labelHashSet keepFaces(2*size());
- forAll(removeFace, facei)
- {
- if (!removeFace[facei])
- {
- keepFaces.insert(facei);
- }
- }
+
+ // Flip from remove face -> keep face
+ faceSelection.flip();
labelList pointMap;
labelList faceMap;
@@ -1425,7 +1445,7 @@ Foam::isoSurfaceTopo::isoSurfaceTopo
(
MeshStorage::subsetMesh
(
- keepFaces,
+ faceSelection,
pointMap,
faceMap
)
diff --git a/src/sampling/surface/isoSurface/isoSurfaceTopo.H b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
index 3e1b8a313593b96815e242496027cac4f1da4536..462dd0a06903105a28d4e387e281fb98cadf350a 100644
--- a/src/sampling/surface/isoSurface/isoSurfaceTopo.H
+++ b/src/sampling/surface/isoSurface/isoSurfaceTopo.H
@@ -99,13 +99,6 @@ class isoSurfaceTopo
void fixTetBasePtIs();
- //- Does any edge of triangle cross iso value?
- bool isTriCut
- (
- const triFace& tri,
- const scalarField& pointValues
- ) const;
-
//- Determine whether cell is cut
cellCutType calcCutType
(