Commit 2689202a authored by Mark Olesen's avatar Mark Olesen Committed by Andrew Heather
Browse files

ENH: support isoSurfaceTopo with cell subsets (#1374)

- treat the faces that would be exposed on a subset as boundary faces
  for the erosion algorithm

STYLE: adjust code for consistency between isoSurfaceCell and isoSurfaceTopo
parent 3eb966d0
......@@ -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);
......
......@@ -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
......
......@@ -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
)
......
......@@ -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
(
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment