isoSurfaceTopo.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2019 OpenFOAM Foundation
    Copyright (C) 2019-2021 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 "isoSurfaceTopo.H"
#include "polyMesh.H"
#include "volFields.H"
#include "edgeHashes.H"
#include "tetCell.H"
#include "tetPointRef.H"
#include "DynamicField.H"
#include "syncTools.H"
#include "uindirectPrimitivePatch.H"
#include "polyMeshTetDecomposition.H"
#include "foamVtkInternalMeshWriter.H"
#include "foamVtkLineWriter.H"
#include "foamVtkSurfaceWriter.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#include "isoSurfaceBaseMethods.H"
defineIsoSurfaceInterpolateMethods(Foam::isoSurfaceTopo);
 
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    defineTypeNameAndDebug(isoSurfaceTopo, 0);
}
 
 
// Get/set snapIndex (0, 1 or 2) at given position
// 0 = no snap
// 1 = snap to first edge end
// 2 = snap to second edge end
// NB: 4 lower bits left free for regular tet-cut information
 
#undef  SNAP_END_VALUE
#undef  SNAP_END_ENCODE
#define SNAP_END_ENCODE(pos, val)   (((val) << (4 + 2 * pos)))
#define SNAP_END_VALUE(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
 
        // 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;
 
        // Add snap index into regular edge cut index
        // Snap to end if less than approx 1% of the distance.
        // - only valid if there is also a corresponding sign change
        #undef  ADD_SNAP_INDEX
        #define ADD_SNAP_INDEX(pos, d1, d2, idx1, idx2)                \
        switch (cutIndex & (idx1 | idx2))                              \
        {                                                              \
            case idx1 : /* first below, second above */                \
            case idx2 : /* first above, second below */                \
                cutIndex |= SNAP_END_ENCODE                            \
                (                                                      \
                    pos,                                               \
                    ((d1 * 100 < d2) ? 1 : (d2 * 100 < d1) ? 2 : 0)    \
                );                                                     \
                break;                                                 \
        }
 
        ADD_SNAP_INDEX(0, p0, p1, 1, 2);    // Edge 0: 0 -> 1
        ADD_SNAP_INDEX(1, p0, p2, 1, 4);    // Edge 1: 0 -> 2
        ADD_SNAP_INDEX(2, p0, p3, 1, 8);    // Edge 2: 0 -> 3
        ADD_SNAP_INDEX(3, p3, p1, 8, 2);    // Edge 3: 3 -> 1
        ADD_SNAP_INDEX(4, p1, p2, 2, 4);    // Edge 4: 1 -> 2
        ADD_SNAP_INDEX(5, p3, p2, 8, 4);    // Edge 5: 3 -> 2
        #undef ADD_SNAP_INDEX
    }
 
    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);
        }
    }
}
 
 
// Return point reference to mesh points or cell-centres
inline static const point& getMeshPointRef
(
    const polyMesh& mesh,
    const label pointi
)
{
    return
    (
        pointi < mesh.nPoints()
      ? mesh.points()[pointi]
      : mesh.cellCentres()[pointi - mesh.nPoints()]
    );
}
 
} // End namespace Foam
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::isoSurfaceTopo::tetCutAddressing::tetCutAddressing
(
    const label nCutCells,
    const bool useSnap,
    const bool useDebugCuts
)
:
    vertsToPointLookup_(12*nCutCells),
    snapVertsLookup_(0),
 
    pointToFace_(10*nCutCells),
    pointFromDiag_(10*nCutCells),
 
    pointToVerts_(10*nCutCells),
    cutPoints_(12*nCutCells),
 
    debugCutTets_(),
    debugCutTetsOn_(useDebugCuts)
{
    // Per cell: 5 pyramids cut, each generating 2 triangles
 
    // Per cell: number of intersected edges:
    //          - four faces cut so 4 mesh edges + 4 face-diagonal edges
    //          - 4 of the pyramid edges
 
    if (useSnap)
    {
        // Some, but not all, cells may have point snapping
        snapVertsLookup_.resize(4*nCutCells);
    }
    if (debugCutTetsOn_)
    {
        debugCutTets_.reserve(6*nCutCells);
    }
}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
void Foam::isoSurfaceTopo::tetCutAddressing::clearDebug()
{
    debugCutTets_.clearStorage();
}
 
 
void Foam::isoSurfaceTopo::tetCutAddressing::clearDiagonal()
{
    pointToFace_.clearStorage();
    pointFromDiag_.clearStorage();
}
 
 
void Foam::isoSurfaceTopo::tetCutAddressing::clearHashes()
{
    vertsToPointLookup_.clear();
    snapVertsLookup_.clear();
}
 
 
Foam::label Foam::isoSurfaceTopo::tetCutAddressing::generatePoint
(
    label facei,
    bool edgeIsDiagonal,
    const int snapEnd,
    const edge& vertices
)
{
    // Generate new point, unless it already exists for edge
    // or corresponds to a snapped point (from another edge)
 
    label pointi = vertsToPointLookup_.lookup(vertices, -1);
    if (pointi == -1)
    {
        bool addNewPoint(true);
 
        const label snapPointi =
        (
            (snapEnd == 1) ? vertices.first()
          : (snapEnd == 2) ? 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;
            edgeIsDiagonal = 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(edgeIsDiagonal);
        }
 
        vertsToPointLookup_.insert(vertices, pointi);
    }
 
    return pointi;
}
 
 
bool Foam::isoSurfaceTopo::tetCutAddressing::generatePoints
(
    const label facei,
    const int tetCutIndex,
    const tetCell& tetLabels,
 
    // Per tet edge whether is face diag etc
    const FixedList<bool, 6>& edgeIsDiagonal
)
{
    bool flip(false);
    const label nCutPointsOld(cutPoints_.size());
 
    // Form the vertices of the triangles for each case
    switch (tetCutIndex & 0x0F)
    {
        case 0x00:
        case 0x0F:
        break;
 
        // Cut point 0
        case 0x0E: flip = true; [[fallthrough]];    // Point 0 above cut
        case 0x01:                                  // Point 0 below cut
        {
            const label cutA
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[0],
                    SNAP_END_VALUE(0, tetCutIndex),
                    tetLabels.edge(0)  // 0 -> 1
                )
            );
            const label cutB
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[1],
                    SNAP_END_VALUE(1, tetCutIndex),
                    tetLabels.edge(1)  // 0 -> 2
                )
            );
            const label cutC
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[2],
                    SNAP_END_VALUE(2, tetCutIndex),
                    tetLabels.edge(2)  // 0 -> 3
                )
            );
 
            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
        }
        break;
 
        // Cut point 1
        case 0x0D: flip = true; [[fallthrough]];    // Point 1 above cut
        case 0x02:                                  // Point 1 below cut
        {
            const label cutA
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[0],
                    SNAP_END_VALUE(0, tetCutIndex),
                    tetLabels.edge(0)  // 0 -> 1
                )
            );
            const label cutB
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[3],
                    SNAP_END_VALUE(3, tetCutIndex),
                    tetLabels.edge(3)  // 3 -> 1
                )
            );
            const label cutC
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[4],
                    SNAP_END_VALUE(4, tetCutIndex),
                    tetLabels.edge(4)  // 1 -> 2
                )
            );
 
            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
        }
        break;
 
        // 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 cutA
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[1],
                    SNAP_END_VALUE(1, tetCutIndex),
                    tetLabels.edge(1)  // 0 -> 2
                )
            );
            const label cutB
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[2],
                    SNAP_END_VALUE(2, tetCutIndex),
                    tetLabels.edge(2)  // 0 -> 3
                )
            );
            const label cutC
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[3],
                    SNAP_END_VALUE(3, tetCutIndex),
                    tetLabels.edge(3)  // 3 -> 1
                )
            );
            const label cutD
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[4],
                    SNAP_END_VALUE(4, tetCutIndex),
                    tetLabels.edge(4)  // 1 -> 2
                )
            );
 
            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
            appendTriLabels(cutPoints_, cutA, cutC, cutD, flip);
        }
        break;
 
        // Cut point 2
        case 0x0B: flip = true; [[fallthrough]];    // Point 2 above cut
        case 0x04:                                  // Point 2 below cut
        {
            const label cutA
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[1],
                    SNAP_END_VALUE(1, tetCutIndex),
                    tetLabels.edge(1)  // 0 -> 2
                )
            );
            const label cutB
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[4],
                    SNAP_END_VALUE(4, tetCutIndex),
                    tetLabels.edge(4)  // 1 -> 2
                )
            );
            const label cutC
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[5],
                    SNAP_END_VALUE(5, tetCutIndex),
                    tetLabels.edge(5)  // 3 -> 2
                )
            );
 
            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
        }
        break;
 
        // 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 cutA
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[0],
                    SNAP_END_VALUE(0, tetCutIndex),
                    tetLabels.edge(0)  // 0 -> 1
                )
            );
            const label cutB
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[4],
                    SNAP_END_VALUE(4, tetCutIndex),
                    tetLabels.edge(4)  // 1 -> 2
                )
            );
            const label cutC
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[5],
                    SNAP_END_VALUE(5, tetCutIndex),
                    tetLabels.edge(5)  // 3 -> 2
                )
            );
            const label cutD
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[2],
                    SNAP_END_VALUE(2, tetCutIndex),
                    tetLabels.edge(2)  // 0 -> 3
                )
            );
 
            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
            appendTriLabels(cutPoints_, cutA, cutC, cutD, flip);
        }
        break;
 
        // 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 cutA
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[0],
                    SNAP_END_VALUE(0, tetCutIndex),
                    tetLabels.edge(0)  // 0 -> 1
                )
            );
            const label cutB
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[3],
                    SNAP_END_VALUE(3, tetCutIndex),
                    tetLabels.edge(3)  // 3 -> 1
                )
            );
            const label cutC
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[5],
                    SNAP_END_VALUE(5, tetCutIndex),
                    tetLabels.edge(5)  // 3 -> 2
                )
            );
            const label cutD
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[1],
                    SNAP_END_VALUE(1, tetCutIndex),
                    tetLabels.edge(1)  // 0 -> 2
                )
            );
 
            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
            appendTriLabels(cutPoints_, cutA, cutC, cutD, flip);
        }
        break;
 
        // Cut point 3
        case 0x07: flip = true; [[fallthrough]];    // Point 3 above cut
        case 0x08:                                  // Point 3 below cut
        {
            const label cutA
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[2],
                    SNAP_END_VALUE(2, tetCutIndex),
                    tetLabels.edge(2)  // 0 -> 3
                )
            );
            const label cutB
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[5],
                    SNAP_END_VALUE(5, tetCutIndex),
                    tetLabels.edge(5)  // 3 -> 2
                )
            );
            const label cutC
            (
                generatePoint
                (
                    facei,
                    edgeIsDiagonal[3],
                    SNAP_END_VALUE(3, tetCutIndex),
                    tetLabels.edge(3)  // 3 -> 1
                )
            );
 
            appendTriLabels(cutPoints_, cutA, cutB, cutC, flip);
        }
        break;
    }
 
    const bool added(nCutPointsOld != cutPoints_.size());
 
    if (added && debugCutTetsOn_)
    {
        debugCutTets_.append(tetLabels.shape());
    }
 
    return added;
}
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
// Requires mesh_, tetBasePtIs
void Foam::isoSurfaceTopo::generateTriPoints
(
    const label celli,
    const bool isTet,
    const labelList& tetBasePtIs,
    tetCutAddressing& tetCutAddr
) 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 facei = cFaces[0];
        const face& f0 = faces[facei];
 
        // Get the other point from f1. Tbd: check if not duplicate face
        // (ACMI / ignoreBoundaryFaces_).
        const face& f1 = faces[cFaces[1]];
        label apexi = -1;
        forAll(f1, fp)
        {
            apexi = f1[fp];
            if (!f0.found(apexi))
            {
                break;
            }
        }
 
        const label p0 = f0[0];
        label p1 = f0[1];
        label p2 = f0[2];
 
        if (faceOwner[facei] == celli)
        {
            std::swap(p1, p2);
        }
 
        const tetCell tetLabels(p0, p1, p2, apexi);
        const int tetCutIndex
        (
            getTetCutIndex
            (
                pVals_[p0],
                pVals_[p1],
                pVals_[p2],
                pVals_[apexi],
                iso_,
                doSnap
            )
        );
 
        tetCutAddr.generatePoints
        (
            facei,
            tetCutIndex,
            tetLabels,
            FixedList<bool, 6>(false)  // Not face diagonal
        );
    }
    else
    {
        for (const label facei : cFaces)
        {
            if
            (
               !mesh_.isInternalFace(facei)
             && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
            )
            {
                continue;
            }
 
            const face& f = faces[facei];
 
            label fp0 = tetBasePtIs[facei];
 
            // Fallback
            if (fp0 < 0)
            {
                fp0 = 0;
            }
 
            const label p0 = f[fp0];
            label fp = f.fcIndex(fp0);
            for (label i = 2; i < f.size(); ++i)
            {
                label p1 = f[fp];
                fp = f.fcIndex(fp);
                label p2 = f[fp];
 
                FixedList<bool, 6> edgeIsDiagonal(false);
                if (faceOwner[facei] == celli)
                {
                    std::swap(p1, p2);
                    if (i != 2) edgeIsDiagonal[1] = true;
                    if (i != f.size()-1) edgeIsDiagonal[0] = true;
                }
                else
                {
                    if (i != 2) edgeIsDiagonal[0] = true;
                    if (i != f.size()-1) edgeIsDiagonal[1] = true;
                }
 
                const tetCell tetLabels(p0, p1, p2, mesh_.nPoints()+celli);
                const int tetCutIndex
                (
                    getTetCutIndex
                    (
                        pVals_[p0],
                        pVals_[p1],
                        pVals_[p2],
                        cVals_[celli],
                        iso_,
                        doSnap
                    )
                );
 
                tetCutAddr.generatePoints
                (
                    facei,
                    tetCutIndex,
                    tetLabels,
                    edgeIsDiagonal
                );
            }
        }
    }
}
 
 
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
void Foam::isoSurfaceTopo::triangulateOutside
(
    const bool filterDiag,
    const primitivePatch& pp,
    const boolUList& pointFromDiag,
    const labelUList& pointToFace,
    const label cellID,
 
    // outputs
    DynamicList<face>& compactFaces,
    DynamicList<label>& compactCellIDs
)
{
    // We can form pockets:
    // - 1. triangle on face
    // - 2. multiple triangles on interior (from diag edges)
    // - the edge loop will be pocket since it is only the diag
    //   edges that give it volume?
 
    // Retriangulate the exterior loops
    const labelListList& edgeLoops = pp.edgeLoops();
    const labelList& mp = pp.meshPoints();
 
    for (const labelList& loop : edgeLoops)
    {
        if (loop.size() > 2)
        {
            compactFaces.append(face(loop.size()));
            face& f = compactFaces.last();
 
            label fpi = 0;
            forAll(f, i)
            {
                const label pointi = mp[loop[i]];
                if (filterDiag && pointFromDiag[pointi])
                {
                    const label prevPointi = mp[loop[loop.fcIndex(i)]];
                    if
                    (
                        pointFromDiag[prevPointi]
                     && (pointToFace[pointi] != pointToFace[prevPointi])
                    )
                    {
                        f[fpi++] = pointi;
                    }
                    else
                    {
                        // Filter out diagonal point
                    }
                }
                else
                {
                    f[fpi++] = pointi;
                }
            }
 
            if (fpi > 2)
            {
                f.resize(fpi);
            }
            else
            {
                // Keep original face
                forAll(f, i)
                {
                    const label pointi = mp[loop[i]];
                    f[i] = pointi;
                }
            }
            compactCellIDs.append(cellID);
        }
    }
}
 
 
void Foam::isoSurfaceTopo::removeInsidePoints
(
    Mesh& s,
    const bool filterDiag,
 
    // inputs
    const boolUList& pointFromDiag,
    const labelUList& pointToFace,
    const labelUList& start,                // Per cell the starting triangle
 
    // outputs
    DynamicList<label>& compactCellIDs      // Per returned tri the cellID
)
{
    const pointField& points = s.points();
 
    compactCellIDs.clear();
    compactCellIDs.reserve(s.size()/4);
 
    DynamicList<face> compactFaces(s.size()/4);
 
    for (label celli = 0; celli < start.size()-1; ++celli)
    {
        // Triangles for the current cell
 
        const label nTris = start[celli+1]-start[celli];
 
        if (nTris)
        {
            const primitivePatch pp
            (
                SubList<face>(s, nTris, start[celli]),
                points
            );
 
            triangulateOutside
            (
                filterDiag,
                pp,
                pointFromDiag,
                pointToFace,
                celli,
 
                compactFaces,
                compactCellIDs
            );
        }
    }
 
    s.swapFaces(compactFaces);  // Use new faces
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::isoSurfaceTopo::isoSurfaceTopo
(
    const polyMesh& mesh,
    const scalarField& cellValues,
    const scalarField& pointValues,
    const scalar iso,
    const isoSurfaceParams& params,
    const bitSet& ignoreCells
)
:
    isoSurfaceBase(mesh, cellValues, pointValues, iso, params)
{
    // The cell cut type
    List<cutType> cellCutType_(mesh.nCells(), cutType::UNVISITED);
 
    // Time description (for debug output)
    const word timeDesc(word::printf("%08d", mesh_.time().timeIndex()));
 
    if (debug)
    {
        Pout<< "isoSurfaceTopo:" << nl
            << "    cell min/max  : " << minMax(cVals_) << nl
            << "    point min/max : " << minMax(pVals_) << nl
            << "    isoValue      : " << iso << nl
            << "    filter        : "
            << isoSurfaceParams::filterNames[params.filter()] << nl
            << "    mesh span     : " << mesh.bounds().mag() << nl
            << "    ignoreCells   : " << ignoreCells.count()
            << " / " << cVals_.size() << nl
            << endl;
    }
 
    this->ignoreCyclics();
 
    label nBlockedCells = 0;
 
    // Mark ignoreCells as blocked
    nBlockedCells += blockCells(cellCutType_, ignoreCells);
 
    // Mark cells outside bounding box as blocked
    nBlockedCells +=
        blockCells(cellCutType_, params.getClipBounds(), volumeType::OUTSIDE);
 
    // Adjusted tet base points to improve tet quality
    labelList tetBasePtIs
    (
        polyMeshTetDecomposition::adjustTetBasePtIs(mesh_, debug)
    );
 
 
    // Determine cell cuts
    const label nCutCells = calcCellCuts(cellCutType_);
 
    if (debug)
    {
        Pout<< "isoSurfaceTopo : candidate cells cut "
            << nCutCells
            << " blocked " << nBlockedCells
            << " total " << mesh_.nCells() << endl;
    }
 
    if (debug && isA<fvMesh>(mesh))
    {
        const auto& fvmesh = dynamicCast<const fvMesh>(mesh);
 
        volScalarField debugField
        (
            IOobject
            (
                "isoSurfaceTopo.cutType",
                fvmesh.time().timeName(),
                fvmesh.time(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            fvmesh,
            dimensionedScalar(dimless)
        );
 
        auto& debugFld = debugField.primitiveFieldRef();
 
        forAll(cellCutType_, celli)
        {
            debugFld[celli] = cellCutType_[celli];
        }
 
        Info<< "Writing cut types: " << debugField.objectRelPath() << nl;
        debugField.write();
    }
 
    // Additional debugging
    if (debug & 8)
    {
        // Write debug cuts cells in VTK format
        {
            constexpr uint8_t realCut(cutType::CUT | cutType::TETCUT);
            labelList debugCutCells(nCutCells, Zero);
 
            label nout = 0;
            forAll(cellCutType_, celli)
            {
                if ((cellCutType_[celli] & realCut) != 0)
                {
                    debugCutCells[nout] = celli;
                    ++nout;
                    if (nout >= nCutCells) break;
                }
            }
 
            // The mesh subset cut
            vtk::vtuCells vtuCells;
            vtuCells.reset(mesh_, debugCutCells);
 
            vtk::internalMeshWriter writer
            (
                mesh_,
                vtuCells,
                fileName
                (
                    mesh_.time().globalPath()
                  / ("isoSurfaceTopo." + timeDesc + "-cutCells")
                )
            );
 
            writer.writeGeometry();
 
            // CellData
            writer.beginCellData();
            writer.writeProcIDs();
            writer.writeCellData("cutField", cVals_);
 
            // PointData
            writer.beginPointData();
            writer.writePointData("cutField", pVals_);
 
            Info<< "isoSurfaceTopo : (debug) wrote "
                << returnReduce(nCutCells, sumOp<label>())
                << " cut cells: "
                << writer.output().name() << nl;
        }
    }
 
 
    tetCutAddressing tetCutAddr
    (
        nCutCells,
        this->snap(),
        (debug & 8)  // Enable debug tets
    );
 
    labelList startTri(mesh_.nCells()+1, Zero);
    for (label celli = 0; celli < mesh_.nCells(); ++celli)
    {
        startTri[celli] = tetCutAddr.nFaces();
        if ((cellCutType_[celli] & cutType::ANYCUT) != 0)
        {
            generateTriPoints
            (
                celli,
                // Same as tetMatcher::test(mesh_, celli),
                bool(cellCutType_[celli] & cutType::TETCUT),
 
                tetBasePtIs,
                tetCutAddr
            );
        }
    }
    startTri.last() = tetCutAddr.nFaces();
 
    // Information not needed anymore:
    tetBasePtIs.clear();
    tetCutAddr.clearHashes();
 
 
    // From list of vertices -> triangular faces
    faceList allTriFaces(startTri.last());
    {
        auto& verts = tetCutAddr.cutPoints();
 
        label verti = 0;
        for (face& f : allTriFaces)
        {
            f.resize(3);
            f[0] = verts[verti++];
            f[1] = verts[verti++];
            f[2] = verts[verti++];
        }
        verts.clearStorage();  // Not needed anymore
    }
 
 
    // The cells cut by the triangular faces
    meshCells_.resize(startTri.last());
    for (label celli = 0; celli < startTri.size()-1; ++celli)
    {
        // All triangles for the current cell
        labelList::subList
        (
            meshCells_,
            (startTri[celli+1] - startTri[celli]),
            startTri[celli]
        ) = celli;
    }
 
 
    pointToVerts_.transfer(tetCutAddr.pointToVerts());
 
    pointField allTriPoints
    (
        this->interpolateTemplate
        (
            mesh_.cellCentres(),
            mesh_.points()
        )
    );
 
 
    // Assign to MeshedSurface
    static_cast<Mesh&>(*this) = Mesh
    (
        std::move(allTriPoints),
        std::move(allTriFaces),
        surfZoneList()  // zones not required (one zone)
    );
 
    if (debug)
    {
        Pout<< "isoSurfaceTopo : generated "
            << Mesh::size() << " triangles "
            << Mesh::points().size() << " points" << endl;
    }
 
    // Write debug triangulated surface
    if ((debug & 8) && (params.filter() != filterType::NONE))
    {
        const Mesh& s = *this;
 
        vtk::surfaceWriter writer
        (
            s.points(),
            s,
            fileName
            (
                mesh_.time().globalPath()
              / ("isoSurfaceTopo." + timeDesc + "-triangles")
            )
        );
 
        writer.writeGeometry();
 
        // CellData
        writer.beginCellData();
        writer.writeProcIDs();
        writer.write("cellID", meshCells_);
 
        // PointData
        writer.beginPointData();
        {
            // NB: may have non-compact surface points
            // --> use points().size() not nPoints()!
 
            labelList pointStatus(s.points().size(), Zero);
 
            forAll(pointToVerts_, i)
            {
                const edge& verts = pointToVerts_[i];
                if (verts.first() == verts.second())
                {
                    // Duplicate index (ie, snapped)
                    pointStatus[i] = 1;
                }
                if (tetCutAddr.pointFromDiag().test(i))
                {
                    // Point on triangulation diagonal
                    pointStatus[i] = -1;
                }
            }
 
            writer.write("point-status", pointStatus);
        }
 
        Info<< "isoSurfaceTopo : (debug) wrote "
            << returnReduce(s.size(), sumOp<label>())
            << " triangles : "
            << writer.output().name() << nl;
    }
 
 
    // Now:
    // - 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
    //    (note that the pyramid faces are shared between multiple mesh faces)
    //  - pointToVerts_ : originating mesh vertex or cell centre
 
    if (params.filter() == filterType::NONE)
    {
        // Compact out unused (snapped) points
        if (this->snap())
        {
            Mesh& s = *this;
 
            labelList pointMap;          // Back to original point
            s.compactPoints(pointMap);   // Compact out unused points
            pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
        }
    }
    else
    {
        // Initial filtering
 
        Mesh& s = *this;
 
        // Triangulate outside
        // (filter edges to cell centres and optionally face diagonals)
        DynamicList<label> compactCellIDs;  // Per tri the cell
 
        removeInsidePoints
        (
            *this,
            // Filter face diagonal
            (
                params.filter() == filterType::DIAGCELL
             || params.filter() == filterType::NONMANIFOLD
            ),
            tetCutAddr.pointFromDiag(),
            tetCutAddr.pointToFace(),
            startTri,
            compactCellIDs
        );
 
        labelList pointMap;          // Back to original point
        s.compactPoints(pointMap);   // Compact out unused points
 
        pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
        meshCells_.transfer(compactCellIDs);
 
        if (debug)
        {
            Pout<< "isoSurfaceTopo :"
                " after removing cell centre and face-diag triangles: "
                << Mesh::size() << " faces "
                << Mesh::points().size() << " points"
                << endl;
        }
    }
 
    // Diagonal filter information not needed anymore
    tetCutAddr.clearDiagonal();
 
 
    // For more advanced filtering (eg, removal of open edges)
    // need the boundary and other 'protected' points
 
    bitSet isProtectedPoint;
    if
    (
        (params.filter() == filterType::NONMANIFOLD)
     || tetCutAddr.debugCutTetsOn()
    )
    {
        // Mark points on mesh outside as 'protected'
        // - never erode these edges
 
        isProtectedPoint.resize(mesh_.nPoints());
 
        for
        (
            label facei = mesh_.nInternalFaces();
            facei < mesh_.nFaces();
            ++facei
        )
        {
            isProtectedPoint.set(mesh_.faces()[facei]);
        }
 
        // Include faces that would be exposed from mesh subset
        if (nBlockedCells)
        {
            const labelList& faceOwn = mesh_.faceOwner();
            const labelList& faceNei = mesh_.faceNeighbour();
 
            for (label facei = 0; facei < mesh.nInternalFaces(); ++facei)
            {
                // If only one cell is blocked, the face corresponds
                // to an exposed subMesh face
                if
                (
                    (cellCutType_[faceOwn[facei]] == cutType::BLOCKED)
                 != (cellCutType_[faceNei[facei]] == cutType::BLOCKED)
                )
                {
                    isProtectedPoint.set(mesh_.faces()[facei]);
                }
            }
        }
    }
 
    // Initial cell cut information not needed anymore
    cellCutType_.clear();
 
 
    // Additional debugging
    if (tetCutAddr.debugCutTetsOn())
    {
        // Write debug cut tets in VTK format
        {
            const auto& debugCuts = tetCutAddr.debugCutTets();
 
            // The TET shapes, using the mesh_ points information
            vtk::vtuCells vtuCells;
            vtuCells.resetShapes(debugCuts);
 
            // Use all points and all cell-centres
            vtuCells.setNumPoints(mesh_.nPoints());
            vtuCells.addPointCellLabels(identity(mesh_.nCells()));
 
            vtk::internalMeshWriter writer
            (
                mesh_,
                vtuCells,
                fileName
                (
                    mesh_.time().globalPath()
                  / ("isoSurfaceTopo." + timeDesc + "-cutTets")
                )
            );
 
            writer.writeGeometry();
 
            // CellData
            writer.beginCellData();
            writer.writeProcIDs();
 
            // Quality of the cut tets
            {
                Field<scalar> cutTetQuality(debugCuts.size());
                forAll(cutTetQuality, teti)
                {
                    cutTetQuality[teti] = tetPointRef
                    (
                        getMeshPointRef(mesh_, debugCuts[teti][0]),
                        getMeshPointRef(mesh_, debugCuts[teti][1]),
                        getMeshPointRef(mesh_, debugCuts[teti][2]),
                        getMeshPointRef(mesh_, debugCuts[teti][3])
                    ).quality();
                }
                writer.writeCellData("tetQuality", cutTetQuality);
            }
 
            // PointData
            if (this->snap())
            {
                writer.beginPointData();
 
                labelList pointStatus(vtuCells.nFieldPoints(), Zero);
 
                for (const edge& verts : pointToVerts_)
                {
                    if (verts.first() == verts.second())
                    {
                        // Duplicate index (ie, snapped)
                        pointStatus[verts.first()] = 1;
                    }
                }
 
                writer.writePointData("point-status", pointStatus);
            }
 
            Info<< "isoSurfaceTopo : (debug) wrote "
                << returnReduce(debugCuts.size(), sumOp<label>())
                << " cut tets: "
                << writer.output().name() << nl;
        }
 
        // Determining open edges. Same logic as used later...
 
        labelHashSet openEdgeIds(0);
 
        {
            const Mesh& s = *this;
 
            const labelList& mp = s.meshPoints();
            const edgeList& surfEdges = s.edges();
            const labelListList& edgeFaces = s.edgeFaces();
            openEdgeIds.resize(2*s.size());
 
            forAll(edgeFaces, edgei)
            {
                const labelList& eFaces = edgeFaces[edgei];
                if (eFaces.size() == 1)
                {
                    // 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()]];
 
                    if
                    (
                        isProtectedPoint.test(verts0.first())
                     && isProtectedPoint.test(verts0.second())
                     && isProtectedPoint.test(verts1.first())
                     && isProtectedPoint.test(verts1.second())
                    )
                    {
                        // Open edge on boundary face. Keep
                    }
                    else
                    {
                        // Open edge
                        openEdgeIds.insert(edgei);
                    }
                }
            }
 
            const label nOpenEdges
            (
                returnReduce(openEdgeIds.size(), sumOp<label>())
            );
 
            if (nOpenEdges)
            {
                const edgeList debugEdges
                (
                    surfEdges,
                    openEdgeIds.sortedToc()
                );
 
                vtk::lineWriter writer
                (
                    s.points(),
                    debugEdges,
                    fileName
                    (
                        mesh_.time().globalPath()
                      / ("isoSurfaceTopo." + timeDesc + "-openEdges")
                    )
                );
 
                writer.writeGeometry();
 
                // CellData
                writer.beginCellData();
                writer.writeProcIDs();
 
                Info<< "isoSurfaceTopo : (debug) wrote "
                    << nOpenEdges << " open edges: "
                    << writer.output().name() << nl;
            }
            else
            {
                Info<< "isoSurfaceTopo : no open edges" << nl;
            }
        }
 
        // Write debug surface with snaps
        if (this->snap())
        {
            const Mesh& s = *this;
 
            vtk::surfaceWriter writer
            (
                s.points(),
                s,
                fileName
                (
                    mesh_.time().globalPath()
                  / ("isoSurfaceTopo." + timeDesc + "-surface")
                )
            );
 
            writer.writeGeometry();
 
            // CellData
            writer.beginCellData();
            writer.writeProcIDs();
            writer.write("cellID", meshCells_);
 
            // PointData
            writer.beginPointData();
            {
                // NB: may have non-compact surface points
                // --> use points().size() not nPoints()!
 
                labelList pointStatus(s.points().size(), Zero);
 
                forAll(pointToVerts_, i)
                {
                    const edge& verts = pointToVerts_[i];
                    if (verts.first() == verts.second())
                    {
                        // Duplicate index (ie, snapped)
                        pointStatus[i] = 1;
                    }
                }
 
                writer.write("point-status", pointStatus);
            }
 
            Info<< "isoSurfaceTopo : (debug) wrote "
                << returnReduce(s.size(), sumOp<label>())
                << " faces : "
                << writer.output().name() << nl;
        }
    }
    tetCutAddr.clearDebug();
 
 
    if (params.filter() == filterType::NONMANIFOLD)
    {
        // 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
        // multiply-connected triangles
 
        // Solved by eroding open-edges
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
        // The list of surface faces that should be retained after erosion
        Mesh& surf = *this;
        labelList faceAddr(identity(surf.size()));
 
        bitSet faceSelection;
 
        while (true)
        {
            // Shadow the surface for the purposes of erosion
            uindirectPrimitivePatch erosion
            (
                UIndirectList<face>(surf, faceAddr),
                surf.points()
            );
 
            faceSelection.clear();
            faceSelection.resize(erosion.size());
 
            const labelList& mp = erosion.meshPoints();
            const edgeList& surfEdges = erosion.edges();
            const labelListList& edgeFaces = erosion.edgeFaces();
 
            label nEdgeRemove = 0;
 
            forAll(edgeFaces, edgei)
            {
                const labelList& eFaces = edgeFaces[edgei];
                if (eFaces.size() == 1)
                {
                    // 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()]];
 
                    if
                    (
                        isProtectedPoint.test(verts0.first())
                     && isProtectedPoint.test(verts0.second())
                     && isProtectedPoint.test(verts1.first())
                     && isProtectedPoint.test(verts1.second())
                    )
                    {
                        // Open edge on boundary face. Keep
                    }
                    else
                    {
                        // Open edge. Mark for erosion
                        faceSelection.set(eFaces[0]);
                        ++nEdgeRemove;
                    }
                }
            }
 
            if (debug)
            {
                Pout<< "isoSurfaceTopo :"
                    << " removing " << faceSelection.count()
                    << " / " << faceSelection.size()
                    << " faces on " << nEdgeRemove << " open edges" << endl;
            }
 
            if (returnReduce(faceSelection.none(), andOp<bool>()))
            {
                break;
            }
 
            // Remove the faces from the addressing
            inplaceSubset(faceSelection, faceAddr, true);  // True = remove
        }
 
 
        // Finished erosion (if any)
        // - retain the faces listed in the updated addressing
 
        if (surf.size() != faceAddr.size())
        {
            faceSelection.clear();
            faceSelection.resize(surf.size());
            faceSelection.set(faceAddr);
 
            inplaceSubsetMesh(faceSelection);
        }
    }
}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
void Foam::isoSurfaceTopo::inplaceSubsetMesh(const bitSet& include)
{
    labelList pointMap;
    labelList faceMap;
    Mesh filtered
    (
        Mesh::subsetMesh(include, pointMap, faceMap)
    );
    Mesh::transfer(filtered);
 
    meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
 
    pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
}
 
 
// ************************************************************************* //