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 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 "isoSurface.H"
#include "polyMesh.H"
#include "tetMatcher.H"
#include "tetPointRef.H"
#include "DynamicField.H"
#include "syncTools.H"
#include "polyMeshTetDecomposition.H"
#include "cyclicACMIPolyPatch.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
    defineTypeNameAndDebug(isoSurfaceTopo, 0);
}
 
 
// * * * * * * * * * * * * * * * 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::isoSurfaceTopo::cellCutType Foam::isoSurfaceTopo::calcCutType
(
    const bool isTet,
    const label celli
) const
{
    if (ignoreCells_.test(celli))
    {
        return NOTCUT;
    }
 
    const cell& cFaces = mesh_.cells()[celli];
 
    if (isTet)
    {
        for (const label facei : cFaces)
        {
            if
            (
               !mesh_.isInternalFace(facei)
             && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
            )
            {
                continue;
            }
 
            const face& f = mesh_.faces()[facei];
 
            for (label fp = 1; fp < f.size() - 1; ++fp)
            {
                if (isTriCut(f[0], f[fp], f[f.fcIndex(fp)], pVals_, iso_))
                {
                    return CUT;
                }
            }
        }
        return NOTCUT;
    }
 
 
    const bool cellLower = (cVals_[celli] < iso_);
 
    // First check if there is any cut in cell
    bool edgeCut = false;
 
    for (const label facei : cFaces)
    {
        if
        (
           !mesh_.isInternalFace(facei)
         && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
        )
        {
            continue;
        }
 
        const face& f = mesh_.faces()[facei];
 
        // Check pyramids cut
        for (const label pointi : f)
        {
            if ((pVals_[pointi] < iso_) != cellLower)
            {
                edgeCut = true;
                break;
            }
        }
 
        if (edgeCut)
        {
            break;
        }
 
        // Check (triangulated) face edges
        // With fallback for problem decompositions
        const label fp0 = (tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei]);
 
        label fp = f.fcIndex(fp0);
        for (label i = 2; i < f.size(); ++i)
        {
            const label nextFp = f.fcIndex(fp);
 
            if (isTriCut(f[fp0], f[fp], f[nextFp], pVals_, iso_))
            {
                edgeCut = true;
                break;
            }
 
            fp = nextFp;
        }
 
        if (edgeCut)
        {
            break;
        }
    }
 
    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
 
        const labelList& cPoints = mesh_.cellPoints(celli);
 
        label nPyrCuts = 0;
 
        for (const label pointi : cPoints)
        {
            if ((pVals_[pointi] < iso_) != cellLower)
            {
                ++nPyrCuts;
            }
        }
 
        if (nPyrCuts == cPoints.size())
        {
            return SPHERE;
        }
        else if (nPyrCuts)
        {
            return CUT;
        }
    }
 
    return NOTCUT;
}
 
 
Foam::label Foam::isoSurfaceTopo::calcCutTypes
(
    tetMatcher& tet,
    List<cellCutType>& cellCutTypes
)
{
    cellCutTypes.setSize(mesh_.nCells());
 
    label nCutCells = 0;
    forAll(cellCutTypes, celli)
    {
        cellCutTypes[celli] = calcCutType(tet.isA(mesh_, celli), celli);
 
        if (cellCutTypes[celli] == CUT)
        {
            ++nCutCells;
        }
    }
 
    if (debug)
    {
        Pout<< "isoSurfaceCell : candidate cut cells "
            << nCutCells << " / " << mesh_.nCells() << endl;
    }
    return nCutCells;
}
 
 
Foam::scalar Foam::isoSurfaceTopo::minTetQ
(
    const label facei,
    const label faceBasePtI
) const
{
    scalar q = polyMeshTetDecomposition::minQuality
    (
        mesh_,
        mesh_.cellCentres()[mesh_.faceOwner()[facei]],
        facei,
        true,
        faceBasePtI
    );
 
    if (mesh_.isInternalFace(facei))
    {
        q = min
        (
            q,
            polyMeshTetDecomposition::minQuality
            (
                mesh_,
                mesh_.cellCentres()[mesh_.faceNeighbour()[facei]],
                facei,
                false,
                faceBasePtI
            )
        );
    }
    return q;
}
 
 
void Foam::isoSurfaceTopo::fixTetBasePtIs()
{
    // Determine points used by two faces on the same cell
    const cellList& cells = mesh_.cells();
    const faceList& faces = mesh_.faces();
    const labelList& faceOwner = mesh_.faceOwner();
    const labelList& faceNeighbour = mesh_.faceNeighbour();
 
 
    // Get face triangulation base point
    tetBasePtIs_ = mesh_.tetBasePtIs();
 
 
    // Pre-filter: mark all cells with illegal base points
    bitSet problemCells(cells.size());
 
    forAll(tetBasePtIs_, facei)
    {
        if (tetBasePtIs_[facei] == -1)
        {
            problemCells.set(faceOwner[facei]);
 
            if (mesh_.isInternalFace(facei))
            {
                problemCells.set(faceNeighbour[facei]);
            }
        }
    }
 
 
    // Mark all points that are shared by just two faces within an adjacent
    // problem cell as problematic
    bitSet problemPoints(mesh_.points().size());
 
    {
        // Number of times a point occurs in a cell.
        // Used to detect dangling vertices (count = 2)
        Map<label> pointCount;
 
        // Analyse problem cells for points shared by two faces only
        for (const label celli : problemCells)
        {
            pointCount.clear();
 
            for (const label facei : cells[celli])
            {
                for (const label pointi : faces[facei])
                {
                    ++pointCount(pointi);
                }
            }
 
            forAllConstIters(pointCount, iter)
            {
                if (iter.val() == 1)
                {
                    FatalErrorInFunction
                        << "point:" << iter.key()
                        << " at:" << mesh_.points()[iter.key()]
                        << " only used by one face" << nl
                        << exit(FatalError);
                }
                else if (iter.val() == 2)
                {
                    problemPoints.set(iter.key());
                }
            }
        }
    }
 
 
    // For all faces which form a part of a problem-cell, check if the base
    // point is adjacent to any problem points. If it is, re-calculate the base
    // point so that it is not.
    label nAdapted = 0;
    forAll(tetBasePtIs_, facei)
    {
        if
        (
            problemCells[faceOwner[facei]]
         || (mesh_.isInternalFace(facei) && problemCells[faceNeighbour[facei]])
        )
        {
            const face& f = faces[facei];
 
            // Check if either of the points adjacent to the base point is a
            // problem point. If not, the existing base point can be retained.
            const label fp0 = tetBasePtIs_[facei] < 0 ? 0 : tetBasePtIs_[facei];
 
            if
            (
                !problemPoints[f[f.rcIndex(fp0)]]
             && !problemPoints[f[f.fcIndex(fp0)]]
            )
            {
                continue;
            }
 
            // A new base point is required. Pick the point that results in the
            // least-worst tet and which is not adjacent to any problem points.
            scalar maxQ = -GREAT;
            label maxFp = -1;
            forAll(f, fp)
            {
                if
                (
                    !problemPoints[f[f.rcIndex(fp)]]
                 && !problemPoints[f[f.fcIndex(fp)]]
                )
                {
                    const scalar q = minTetQ(facei, fp);
                    if (q > maxQ)
                    {
                        maxQ = q;
                        maxFp = fp;
                    }
                }
            }
 
            if (maxFp != -1)
            {
                // Success! Set the new base point
                tetBasePtIs_[facei] = maxFp;
            }
            else
            {
                // No point was found on face that would not result in some
                // duplicate triangle. Do what? Continue and hope? Spit an
                // error? Silently or noisily reduce the filtering level?
 
                tetBasePtIs_[facei] = 0;
            }
 
            ++nAdapted;
        }
    }
 
 
    if (debug)
    {
        Pout<< "isoSurface : adapted starting point of triangulation on "
            << nAdapted << " faces." << endl;
    }
 
    syncTools::syncFaceList(mesh_, tetBasePtIs_, maxEqOp<label>());
}
 
 
Foam::label Foam::isoSurfaceTopo::generatePoint
(
    const label facei,
    const bool edgeIsDiag,
    const edge& vertices,
 
    DynamicList<edge>& pointToVerts,
    DynamicList<label>& pointToFace,
    DynamicList<bool>& pointFromDiag,
    EdgeMap<label>& vertsToPoint
) const
{
    const auto edgeFnd = vertsToPoint.cfind(vertices);
    if (edgeFnd.found())
    {
        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;
}
 
 
void Foam::isoSurfaceTopo::generateTriPoints
(
    const label facei,
    const FixedList<scalar, 4>& s,
    const FixedList<point, 4>& p,
    const FixedList<label, 4>& pIndex,
    const FixedList<bool, 6>& edgeIsDiag,// Per tet edge whether is face diag
 
    DynamicList<edge>& pointToVerts,
    DynamicList<label>& pointToFace,
    DynamicList<bool>& pointFromDiag,
 
    EdgeMap<label>& vertsToPoint,
    DynamicList<label>& verts       // Every three verts is new triangle
) const
{
    int triIndex = 0;
    if (s[0] < iso_)
    {
        triIndex |= 1;
    }
    if (s[1] < iso_)
    {
        triIndex |= 2;
    }
    if (s[2] < iso_)
    {
        triIndex |= 4;
    }
    if (s[3] < iso_)
    {
        triIndex |= 8;
    }
 
    // Form the vertices of the triangles for each case
    switch (triIndex)
    {
        case 0x00:
        case 0x0F:
        break;
 
        case 0x01:
        case 0x0E:
        {
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[0],
                    edge(pIndex[0], pIndex[1]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[1],
                    edge(pIndex[0], pIndex[2]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[2],
                    edge(pIndex[0], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
 
            if (triIndex == 0x0E)
            {
                // Flip normals
                label sz = verts.size();
                Swap(verts[sz-2], verts[sz-1]);
            }
        }
        break;
 
        case 0x02:
        case 0x0D:
        {
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[0],
                    edge(pIndex[1], pIndex[0]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[3],
                    edge(pIndex[1], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[4],
                    edge(pIndex[1], pIndex[2]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
 
            if (triIndex == 0x0D)
            {
                // Flip normals
                label sz = verts.size();
                Swap(verts[sz-2], verts[sz-1]);
            }
        }
        break;
 
        case 0x03:
        case 0x0C:
        {
            label p0p2
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[1],
                    edge(pIndex[0], pIndex[2]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            label p1p3
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[3],
                    edge(pIndex[1], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
 
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[2],
                    edge(pIndex[0], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append(p1p3);
            verts.append(p0p2);
            verts.append(p1p3);
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[4],
                    edge(pIndex[1], pIndex[2]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append(p0p2);
 
            if (triIndex == 0x0C)
            {
                // Flip normals
                label sz = verts.size();
                Swap(verts[sz-5], verts[sz-4]);
                Swap(verts[sz-2], verts[sz-1]);
            }
        }
        break;
 
        case 0x04:
        case 0x0B:
        {
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[1],
                    edge(pIndex[2], pIndex[0]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[4],
                    edge(pIndex[2], pIndex[1]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[5],
                    edge(pIndex[2], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
 
            if (triIndex == 0x0B)
            {
                // Flip normals
                label sz = verts.size();
                Swap(verts[sz-2], verts[sz-1]);
            }
        }
        break;
 
        case 0x05:
        case 0x0A:
        {
            label p0p1
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[0],
                    edge(pIndex[0], pIndex[1]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            label p2p3
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[5],
                    edge(pIndex[2], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
 
            verts.append(p0p1);
            verts.append(p2p3);
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[2],
                    edge(pIndex[0], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append(p0p1);
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[4],
                    edge(pIndex[1], pIndex[2]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append(p2p3);
 
            if (triIndex == 0x0A)
            {
                // Flip normals
                label sz = verts.size();
                Swap(verts[sz-5], verts[sz-4]);
                Swap(verts[sz-2], verts[sz-1]);
            }
        }
        break;
 
        case 0x06:
        case 0x09:
        {
            label p0p1
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[0],
                    edge(pIndex[0], pIndex[1]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            label p2p3
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[5],
                    edge(pIndex[2], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
 
            verts.append(p0p1);
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[3],
                    edge(pIndex[1], pIndex[3]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append(p2p3);
            verts.append(p0p1);
            verts.append(p2p3);
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[1],
                    edge(pIndex[0], pIndex[2]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
 
            if (triIndex == 0x09)
            {
                // Flip normals
                label sz = verts.size();
                Swap(verts[sz-5], verts[sz-4]);
                Swap(verts[sz-2], verts[sz-1]);
            }
        }
        break;
 
        case 0x08:
        case 0x07:
        {
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[2],
                    edge(pIndex[3], pIndex[0]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[5],
                    edge(pIndex[3], pIndex[2]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            verts.append
            (
                generatePoint
                (
                    facei,
                    edgeIsDiag[3],
                    edge(pIndex[3], pIndex[1]),
                    pointToVerts, pointToFace, pointFromDiag, vertsToPoint
                )
            );
            if (triIndex == 0x07)
            {
                // Flip normals
                label sz = verts.size();
                Swap(verts[sz-2], verts[sz-1]);
            }
        }
        break;
    }
}
 
 
void Foam::isoSurfaceTopo::generateTriPoints
(
    const label celli,
    const bool isTet,
 
    DynamicList<edge>& pointToVerts,
    DynamicList<label>& pointToFace,
    DynamicList<bool>& pointFromDiag,
 
    EdgeMap<label>& vertsToPoint,
    DynamicList<label>& verts,
    DynamicList<label>& faceLabels
) const
{
    const faceList& faces = mesh_.faces();
    const labelList& faceOwner = mesh_.faceOwner();
    const pointField& points = mesh_.points();
    const cell& cFaces = mesh_.cells()[celli];
 
    if (isTet)
    {
        // For tets don't do cell-centre decomposition, just use the
        // tet points and values
 
        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 oppositeI = -1;
        forAll(f1, fp)
        {
            oppositeI = f1[fp];
            if (!f0.found(oppositeI))
            {
                break;
            }
        }
 
 
        label p0 = f0[0];
        label p1 = f0[1];
        label p2 = f0[2];
        FixedList<bool, 6> edgeIsDiag(false);
 
        if (faceOwner[facei] == celli)
        {
            Swap(p1, p2);
        }
 
        tetPointRef tet
        (
            points[p0],
            points[p1],
            points[p2],
            points[oppositeI]
        );
 
        label startTrii = verts.size();
        generateTriPoints
        (
            facei,
            FixedList<scalar, 4>
            ({
                pVals_[p0],
                pVals_[p1],
                pVals_[p2],
                pVals_[oppositeI]
            }),
            FixedList<point, 4>
            ({
                points[p0],
                points[p1],
                points[p2],
                points[oppositeI]
            }),
            FixedList<label, 4>
            ({
                p0,
                p1,
                p2,
                oppositeI
            }),
            edgeIsDiag,
 
            pointToVerts,
            pointToFace,
            pointFromDiag,
            vertsToPoint,
            verts       // Every three verts is new triangle
        );
 
        label nTris = (verts.size()-startTrii)/3;
        for (label i = 0; i < nTris; ++i)
        {
            faceLabels.append(facei);
        }
    }
    else
    {
        const pointField& cellCentres = mesh_.cellCentres();
 
        for (const label facei : cFaces)
        {
            if
            (
               !mesh_.isInternalFace(facei)
             && ignoreBoundaryFaces_.test(facei-mesh_.nInternalFaces())
            )
            {
                continue;
            }
 
            const face& f = faces[facei];
 
            label fp0 = tetBasePtIs_[facei];
 
            label startTrii = verts.size();
 
            // Fallback
            if (fp0 < 0)
            {
                fp0 = 0;
            }
 
            label fp = f.fcIndex(fp0);
            for (label i = 2; i < f.size(); ++i)
            {
                label nextFp = f.fcIndex(fp);
 
                FixedList<bool, 6> edgeIsDiag(false);
 
                label p0 = f[fp0];
                label p1 = f[fp];
                label p2 = f[nextFp];
                if (faceOwner[facei] == celli)
                {
                    Swap(p1, p2);
                    if (i != 2) edgeIsDiag[1] = true;
                    if (i != f.size()-1) edgeIsDiag[0] = true;
                }
                else
                {
                    if (i != 2) edgeIsDiag[0] = true;
                    if (i != f.size()-1) edgeIsDiag[1] = true;
                }
 
                tetPointRef tet
                (
                    points[p0],
                    points[p1],
                    points[p2],
                    cellCentres[celli]
                );
 
                generateTriPoints
                (
                    facei,
                    FixedList<scalar, 4>
                    ({
                        pVals_[p0],
                        pVals_[p1],
                        pVals_[p2],
                        cVals_[celli]
                    }),
                    FixedList<point, 4>
                    ({
                        points[p0],
                        points[p1],
                        points[p2],
                        cellCentres[celli]
                    }),
                    FixedList<label, 4>
                    ({
                        p0,
                        p1,
                        p2,
                        mesh_.nPoints()+celli
                    }),
                    edgeIsDiag,
 
                    pointToVerts,
                    pointToFace,
                    pointFromDiag,
                    vertsToPoint,
                    verts       // Every three verts is new triangle
                );
 
                fp = nextFp;
            }
 
            label nTris = (verts.size()-startTrii)/3;
            for (label i = 0; i < nTris; ++i)
            {
                faceLabels.append(facei);
            }
        }
    }
}
 
 
void Foam::isoSurfaceTopo::triangulateOutside
(
    const bool filterDiag,
    const primitivePatch& pp,
    const boolList& pointFromDiag,
    const labelList& pointToFace,
    const label cellID,
 
    DynamicList<face>& compactFaces,
    DynamicList<label>& compactCellIDs
) const
{
    // 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(0));
            face& f = compactFaces.last();
 
            f.setSize(loop.size());
            label fpi = 0;
            forAll(f, i)
            {
                label pointi = mp[loop[i]];
                if (filterDiag && pointFromDiag[pointi])
                {
                    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.setSize(fpi);
            }
            else
            {
                // Keep original face
                forAll(f, i)
                {
                    label pointi = mp[loop[i]];
                    f[i] = pointi;
                }
            }
            compactCellIDs.append(cellID);
        }
    }
}
 
 
Foam::MeshedSurface<Foam::face> Foam::isoSurfaceTopo::removeInsidePoints
(
    const bool filterDiag,
    const MeshStorage& s,
    const boolList& pointFromDiag,
    const labelList& pointToFace,
    const labelList& start,                 // Per cell the starting triangle
    DynamicList<label>& pointCompactMap,    // Per returned point the original
    DynamicList<label>& compactCellIDs      // Per returned tri the cellID
) const
{
    const pointField& points = s.points();
 
    pointCompactMap.clear();
    compactCellIDs.clear();
 
    DynamicList<face> compactFaces(s.size()/8);
 
    for (label celli = 0; celli < start.size()-1; ++celli)
    {
        // All triangles of the current cell
 
        label nTris = start[celli+1]-start[celli];
 
        if (nTris)
        {
            const SubList<face> cellFaces(s, nTris, start[celli]);
            const primitivePatch pp(cellFaces, points);
 
            triangulateOutside
            (
                filterDiag,
                pp,
                pointFromDiag,
                pointToFace,
                //protectedFace,
                celli,
 
                compactFaces,
                compactCellIDs
            );
        }
    }
 
 
    // Compact out unused points
    // Pick up the used vertices
    labelList oldToCompact(points.size(), -1);
    DynamicField<point> compactPoints(points.size());
    pointCompactMap.clear();
 
    for (face& f : compactFaces)
    {
        forAll(f, fp)
        {
            label pointi = f[fp];
            label compacti = oldToCompact[pointi];
            if (compacti == -1)
            {
                compacti = compactPoints.size();
                oldToCompact[pointi] = compacti;
                compactPoints.append(points[pointi]);
                pointCompactMap.append(pointi);
            }
            f[fp] = compacti;
        }
    }
 
 
    MeshStorage cpSurface
    (
        std::move(compactPoints),
        std::move(compactFaces),
        s.surfZones()
    );
 
    return cpSurface;
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::isoSurfaceTopo::isoSurfaceTopo
(
    const polyMesh& mesh,
    const scalarField& cVals,
    const scalarField& pVals,
    const scalar iso,
    const filterType filter,
    const boundBox& bounds,
    const bitSet& ignoreCells
)
:
    isoSurfaceBase(iso, bounds),
    mesh_(mesh),
    cVals_(cVals),
    pVals_(pVals),
    ignoreCells_(ignoreCells)
{
    if (debug)
    {
        Pout<< "isoSurfaceTopo::"
            << "    cell min/max  : "
            << min(cVals_) << " / "
            << max(cVals_) << nl
            << "    point min/max : "
            << min(pVals_) << " / "
            << max(pVals_) << nl
            << "    isoValue      : " << iso << nl
            << "    filter        : " << isoSurfaceTopo::filterNames[filter]
            << nl
            << "    mesh span     : " << mesh.bounds().mag() << nl
            << "    ignoreCells   : " << ignoreCells_.count()
            << " / " << cVals_.size() << nl
            << endl;
    }
 
    // Determine boundary pyramids to ignore (since originating from ACMI faces)
    // Maybe needs to be optional argument or more general detect colocated
    // faces.
    {
        const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
        forAll(pbm, patchi)
        {
            const polyPatch& pp = pbm[patchi];
            if (isA<cyclicACMIPolyPatch>(pp))
            {
                ignoreBoundaryFaces_.setSize(mesh_.nBoundaryFaces());
                forAll(pp, i)
                {
                    label facei = pp.start()+i;
                    ignoreBoundaryFaces_.set(facei-mesh_.nInternalFaces());
                }
            }
        }
    }
 
 
    fixTetBasePtIs();
 
    tetMatcher tet;
 
    // Determine if any cut through cell
    List<cellCutType> cellCutTypes;
    const label nCutCells = calcCutTypes(tet, cellCutTypes);
 
    // Per cell: 5 pyramids cut, each generating 2 triangles
    //  - pointToVerts : from generated iso point to originating mesh verts
    DynamicList<edge> pointToVerts(10*nCutCells);
    //  - pointToFace : from generated iso point to originating mesh face
    DynamicList<label> pointToFace(10*nCutCells);
    //  - pointToFace : from generated iso point whether 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);
 
 
    labelList startTri(mesh_.nCells()+1, Zero);
 
    for (label celli = 0; celli < mesh_.nCells(); ++celli)
    {
        startTri[celli] = faceLabels.size();
        if (cellCutTypes[celli] != NOTCUT)
        {
            generateTriPoints
            (
                celli,
                tet.isA(mesh_, celli),
 
                pointToVerts,
                pointToFace,
                pointFromDiag,
 
                vertsToPoint,
                verts,
                faceLabels
            );
 
            for (label i = startTri[celli]; i < faceLabels.size(); ++i)
            {
                cellLabels.append(celli);
            }
        }
    }
    startTri[mesh_.nCells()] = faceLabels.size();
 
 
    pointToVerts_.transfer(pointToVerts);
    meshCells_.transfer(cellLabels);
    pointToFace_.transfer(pointToFace);
 
    tmp<pointField> allPoints
    (
        interpolate
        (
            mesh_.cellCentres(),
            mesh_.points()
        )
    );
 
 
    // Assign to MeshedSurface
    faceList allTris(faceLabels.size());
    label verti = 0;
    for (face& allTri : allTris)
    {
        allTri.setSize(3);
        allTri[0] = verts[verti++];
        allTri[1] = verts[verti++];
        allTri[2] = verts[verti++];
    }
 
 
    surfZoneList allZones(one(), surfZone("allFaces", allTris.size()));
 
    MeshStorage::clear();
    MeshStorage updated
    (
        std::move(allPoints),
        std::move(allTris),
        std::move(allZones)
    );
    MeshStorage::transfer(updated);
 
    // 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 (debug)
    {
        Pout<< "isoSurfaceTopo : generated " << size() << " faces." << endl;
    }
 
 
    if (filter != filterType::NONE)
    {
        // Triangulate outside (filter edges to cell centres and optionally
        // face diagonals)
        DynamicList<label> pointCompactMap(size()); // Back to original point
        DynamicList<label> compactCellIDs(size());  // Per tri the cell
        MeshStorage::operator=
        (
            removeInsidePoints
            (
                (filter == filterType::DIAGCELL ? true : false),
                *this,
                pointFromDiag,
                pointToFace_,
                startTri,
                pointCompactMap,
                compactCellIDs
            )
        );
 
        pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointCompactMap)();
        pointToFace_ = UIndirectList<label>(pointToFace_, pointCompactMap)();
        pointFromDiag = UIndirectList<bool>(pointFromDiag, pointCompactMap)();
        meshCells_.transfer(compactCellIDs);
 
        if (debug)
        {
            Pout<< "isoSurfaceTopo :"
                << " after removing cell centre and face-diag triangles : "
                << size() << endl;
        }
 
 
        if (filter == filterType::DIAGCELL)
        {
            // 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
            // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
            // Mark points on mesh outside. Note that we extend with nCells
            // so we can easily index with pointToVerts_.
            bitSet isBoundaryPoint(mesh.nPoints() + mesh.nCells());
            for
            (
                label facei = mesh.nInternalFaces();
                facei < mesh.nFaces();
                ++facei
            )
            {
                isBoundaryPoint.set(mesh.faces()[facei]);
            }
 
 
            // 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();
 
                const labelListList& edgeFaces = MeshStorage::edgeFaces();
 
                forAll(edgeFaces, edgei)
                {
                    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 on boundary face. Keep
                        }
                        else
                        {
                            // Open edge. Mark for erosion
                            faceSelection.set(eFaces[0]);
                        }
                    }
                }
 
                if (debug)
                {
                    Pout<< "isoSurfaceTopo :"
                        << " removing " << faceSelection.count()
                        << " faces since on open edges" << endl;
                }
 
                if (returnReduce(faceSelection.none(), andOp<bool>()))
                {
                    break;
                }
 
 
                // Flip from remove face -> keep face
                faceSelection.flip();
 
                labelList pointMap;
                labelList faceMap;
                MeshStorage filteredSurf
                (
                    MeshStorage::subsetMesh
                    (
                        faceSelection,
                        pointMap,
                        faceMap
                    )
                );
                MeshStorage::transfer(filteredSurf);
 
                pointToVerts_ = UIndirectList<edge>(pointToVerts_, pointMap)();
                pointToFace_ = UIndirectList<label>(pointToFace_, pointMap)();
                pointFromDiag = UIndirectList<bool>(pointFromDiag, pointMap)();
                meshCells_ = UIndirectList<label>(meshCells_, faceMap)();
            }
        }
    }
}
 
 
// ************************************************************************* //