openfoam/src/sampling/surface/cuttingPlane/cuttingPlane.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  | Copyright (C) 2018 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 "cuttingPlane.H"
#include "fvMesh.H"
#include "volFields.H"
#include "meshTools.H"
#include "edgeHashes.H"
#include "HashOps.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

int Foam::cuttingPlane::debug(Foam::debug::debugSwitch("cuttingPlane", 0));


// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //

namespace Foam
{
    // Check edge/plane intersection based on crossings ... trivial check.
    // Orients the edge (first,last) points in the positive normal direction
    inline bool intersectEdgeOrient(const PackedList<2>& sides, edge& e)
    {
        if (sides[e.first()] == sides[e.last()])
        {
            return false;
        }
        else if (sides[e.last()] < sides[e.first()])
        {
            e.flip();
        }

        return true;
    }


    // Classify sides of plane (0=BACK, 1=ONPLANE, 2=FRONT) for each point
    inline PackedList<2> classifySides(const plane& pln, const pointField& pts)
    {
        const label len = pts.size();

        PackedList<2> output(len);

        // From signed (-1,0,+1) to (0,1,2) for PackedList
        for (label i=0; i < len; ++i)
        {
            output.set(i, unsigned(1 + pln.sign(pts[i], SMALL)));
        }

        return output;
    }


    // Check for face/plane intersection based on crossings
    // Took (-1,0,+1) from plane::sign and packed as (0,1,2).
    // Now use for left shift to obtain (1,2,4).
    //
    // Test accumulated value for an intersection with the plane.
    inline bool intersectsFace
    (
        const PackedList<2>& sides,
        const labelUList& indices
    )
    {
        unsigned accum = 0u;

        for (const label pointi : indices)
        {
            accum |= (1u << sides[pointi]);
        }

        // Accumulated value 3,5,6,7 indicates an intersection
        return (accum == 3 || accum >= 5);
    }

} // End namespace Foam


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

Foam::label Foam::cuttingPlane::calcCellCuts
(
    const primitiveMesh& mesh,
    const PackedList<2>& sides,
    bitSet& cellCuts
)
{
    const faceList& faces = mesh.faces();

    const label nCells = mesh.nCells();
    const label nFaces = mesh.nFaces();
    const label nInternFaces = mesh.nInternalFaces();

    // Detect cells cuts from the face cuts

    label nFaceCuts = 0;

    // 1st face cut of cell
    bitSet hasCut1(nCells);

    // 2nd face cut of cell
    bitSet hasCut2(nCells);

    for (label facei = 0; facei < nInternFaces; ++facei)
    {
        if (intersectsFace(sides, faces[facei]))
        {
            const label own = mesh.faceOwner()[facei];
            const label nei = mesh.faceNeighbour()[facei];

            ++nFaceCuts;

            if (!hasCut1.set(own))
            {
                hasCut2.set(own);
            }
            if (!hasCut1.set(nei))
            {
                hasCut2.set(nei);
            }
        }
    }

    for (label facei = nInternFaces; facei < nFaces; ++facei)
    {
        if (intersectsFace(sides, faces[facei]))
        {
            const label own = mesh.faceOwner()[facei];

            ++nFaceCuts;

            if (!hasCut1.set(own))
            {
                hasCut2.set(own);
            }
        }
    }

    hasCut1.clearStorage();   // Not needed now

    if (cellCuts.size())
    {
        cellCuts.resize(nCells);    // safety
        cellCuts &= hasCut2;        // restrict to cell subset

        if (debug)
        {
            Pout<<"detected " << cellCuts.count() << "/" << nCells
                << " cells cut, subsetted from "
                << hasCut2.count() << "/" << nCells << " cells." << endl;
        }
    }
    else
    {
        cellCuts = std::move(hasCut2);

        if (debug)
        {
            Pout<<"detected " << cellCuts.count() << "/" << nCells
                << " cells cut." << endl;
        }
    }


    if (debug && isA<fvMesh>(mesh))
    {
        const fvMesh& fvm = dynamicCast<const fvMesh&>(mesh);

        volScalarField debugField
        (
            IOobject
            (
                "cuttingPlane.cellCuts",
                fvm.time().timeName(),
                fvm.time(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            fvm,
            dimensionedScalar(dimless, Zero)
        );

        auto& debugFld = debugField.primitiveFieldRef();

        for (const label celli : cellCuts)
        {
            debugFld[celli] = 1;
        }

        Pout<< "Writing cut types:"
            << debugField.objectPath() << endl;

        debugField.write();
    }


    return nFaceCuts;
}


void Foam::cuttingPlane::walkCellCuts
(
    const primitiveMesh& mesh,
    const bitSet& cellCuts,
    const PackedList<2>& sides,
    const bool triangulate,
    label nFaceCuts
)
{
    // Information required from the mesh
    const faceList& faces = mesh.faces();
    const cellList& cells = mesh.cells();
    const pointField& points = mesh.points();


    // Dynamic lists to handle triangulation and/or missed cuts etc
    const label nCellCuts = cellCuts.count();

    DynamicList<point> dynCutPoints(4*nCellCuts);
    DynamicList<face>  dynCutFaces(4*nCellCuts);
    DynamicList<label> dynCutCells(nCellCuts);

    // No nFaceCuts provided? Use a reasonable estimate
    if (!nFaceCuts)
    {
        nFaceCuts = 4*nCellCuts;
    }

    // Edge to pointId mapping
    EdgeMap<label> handledEdges(4*nFaceCuts);

    // Local scratch space for face vertices
    DynamicList<label> localFaceLoop(16);

    // Local scratch space for edge to pointId
    EdgeMap<label> localEdges(128);

    // Local scratch space for edge to faceId
    EdgeMap<edge>  localFaces(128);

    // Avoid duplicates for cuts exactly through a mesh point.
    // No other way to distinguish them, since there is no single edge
    // that "owns" the point.
    Map<label> endPoints;

    // Hash of faces (face points) that are exactly on a cell face
    HashSet<labelList, labelList::Hash<>> onCellFace;


    // Failure handling

    // Cells where walking failed (concave, degenerate, ...)
    labelHashSet failedCells;

    // To unwind insertion of end-point cutting
    labelHashSet localEndPoints;

    // Our recovery point on failure
    label unwindPoint = 0;

    // Cleanup routine for failed cell cut:
    const auto unwindWalk =
        [&](const label failedCellId = -1) -> void
        {
            // Discard points introduced
            dynCutPoints.resize(unwindPoint);

            // Discard end-point cuts
            endPoints.erase(localEndPoints);

            // Optionally record the failedCellId
            if (failedCellId != -1)
            {
                failedCells.insert(failedCellId);
            }
        };


    // Loop over cells that are cut
    for (const label celli : cellCuts)
    {
        const cell& cFace = cells[celli];

        // Reset local scratch
        localEdges.clear();
        localFaces.clear();
        localFaceLoop.clear();
        localEndPoints.clear();

        unwindPoint = dynCutPoints.size();


        // Classification for all the points cut - see intersectsFace() above
        // for more detail
        unsigned pointCutType = 0u;

        for (const label facei : cFace)
        {
            const face& f = faces[facei];

            forAll(f, fp)
            {
                edge e(f.faceEdge(fp));

                // Action #1: detect edge intersection and orient edge
                if (!intersectEdgeOrient(sides, e))
                {
                    continue;
                }

                // Record the edge/face combination for the edge cut.
                // NB: the second operation is edge::insert() which places
                // facei in the unoccupied 'slot'
                localFaces(e).insert(facei);

                // Already handled cut point in this inner-loop?
                if (localEdges.found(e))
                {
                    // No new edge cut required
                    continue;
                }

                // Already handled cut point in the outer-loop?
                label cutPointId = handledEdges.lookup(e, -1);

                if (cutPointId >= 0)
                {
                    // Copy existing edge cut-point index
                    localEdges.insert(e, cutPointId);
                    continue;
                }

                // Expected id for the cut point
                cutPointId = dynCutPoints.size();

                const point& p0 = points[e[0]];
                const point& p1 = points[e[1]];

                // Action #2: edge cut alpha
                const scalar alpha =
                    this->lineIntersect(linePointRef(p0, p1));

                if (alpha < SMALL)
                {
                    // -> equal(alpha, 0) with more tolerance

                    if (endPoints.insert(e[0], cutPointId))
                    {
                        localEndPoints.insert(e[0]);
                        dynCutPoints.append(p0);
                    }
                    else
                    {
                        cutPointId = endPoints[e[0]];
                    }

                    pointCutType |= 0x1; // Cut at 0
                }
                else if (alpha >= (1.0 - SMALL))
                {
                    // -> equal(alpha, 1) with more tolerance

                    if (endPoints.insert(e[1], cutPointId))
                    {
                        localEndPoints.insert(e[1]);
                        dynCutPoints.append(p1);
                    }
                    else
                    {
                        cutPointId = endPoints[e[1]];
                    }

                    pointCutType |= 0x2; // Cut at 1
                }
                else
                {
                    dynCutPoints.append((1-alpha)*p0 + alpha*p1);

                    pointCutType |= 0x4; // Cut between
                }

                // Introduce new edge cut point
                localEdges.insert(e, cutPointId);
            }
        }


        // The keys of localEdges, localFaces are now identical.
        // The size of either should represent the number of points for
        // the resulting face loop.

        const label nTargetLoop = localFaces.size();

        if (nTargetLoop < 3)
        {
            unwindWalk(celli);
            continue;
        }


        // Handling cuts between two cells
        // After the previous intersectEdgeOrient call, the edge oriented
        // according to the gradient.
        // If we only ever cut at the same edge end we know that we have
        // a cut coinciding with a cell face.

        if (pointCutType == 1 || pointCutType == 2)
        {
            // Hash the face-points to avoid duplicate faces
            if (!onCellFace.insert(HashTableOps::values(localEdges, true)))
            {
                DebugInfo
                    <<"skip duplicate on-place cut for cell " << celli
                    << " type (" << pointCutType << ")" << endl;

                // A duplicate is not failure
                unwindWalk();
                continue;
            }
        }


        // Start somewhere.

        // Since the local edges are oriented according to the gradient,
        // they can also be used to determine the correct face orientation.

        const edge refEdge = localFaces.begin().key();
        label nextFace = localFaces.begin().object()[0];

        DebugInfo
            << "search face " <<  nextFace << " IN " <<  localEdges << endl;

        for
        (
            label loopi = 0;
            localFaces.size() && loopi < 2*nTargetLoop;
            ++loopi
        )
        {
            bool ok = false;

            forAllIters(localFaces, iter)
            {
                DebugInfo
                    << "lookup " << nextFace << " in " << iter.object() << nl;

                // Find local index (0,1) or -1 on failure
                const label got = iter.object().which(nextFace);

                if (got != -1)
                {
                    ok = true;

                    // The other face
                    nextFace = iter.object()[(got?0:1)];

                    // The edge -> cut point
                    localFaceLoop.append(localEdges[iter.key()]);

                    DebugInfo
                        <<" faces " << iter.object()
                        << " point " << localFaceLoop.last()
                        << " edge=" << iter.key() << " nextFace=" << nextFace
                        << nl;

                    // Done this connection
                    localFaces.erase(iter);
                    break;
                }
            }

            if (!ok)
            {
                break;
            }
        }


        // Could also check if localFaces is now empty.

        if (nTargetLoop != localFaceLoop.size())
        {
            DebugInfo
                << "Warn expected " << nTargetLoop << " but got "
                << localFaceLoop.size() << endl;

            unwindWalk(celli);
            continue;
        }


        // Success
        handledEdges += localEdges;

        face f(localFaceLoop);

        // Orient face to point in the same direction as the edge gradient
        if ((f.areaNormal(dynCutPoints) & refEdge.vec(points)) < 0)
        {
            f.flip();
        }

        // The cut faces can be quite ugly, so optionally triangulate
        if (triangulate)
        {
            label nTri = f.triangles(dynCutPoints, dynCutFaces);
            while (nTri--)
            {
                dynCutCells.append(celli);
            }
        }
        else
        {
            dynCutFaces.append(f);
            dynCutCells.append(celli);
        }
    }

    if (failedCells.size())
    {
        WarningInFunction
            << "Failed cuts for " << failedCells.size() << " cells:" << nl
            << "    " << flatOutput(failedCells.sortedToc()) << nl
            << endl;
    }


    // No cuts? Then no need for any of this information
    if (dynCutCells.empty())
    {
        this->storedPoints().clear();
        this->storedFaces().clear();
        meshCells_.clear();
    }
    else
    {
        this->storedPoints().transfer(dynCutPoints);
        this->storedFaces().transfer(dynCutFaces);
        meshCells_.transfer(dynCutCells);
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::cuttingPlane::cuttingPlane(const plane& pln)
:
    plane(pln)
{}


Foam::cuttingPlane::cuttingPlane
(
    const plane& pln,
    const primitiveMesh& mesh,
    const bool triangulate,
    const bitSet& cellIdLabels
)
:
    plane(pln)
{
    performCut(mesh, triangulate, cellIdLabels);
}


Foam::cuttingPlane::cuttingPlane
(
    const plane& pln,
    const primitiveMesh& mesh,
    const bool triangulate,
    bitSet&& cellIdLabels
)
:
    plane(pln)
{
    performCut(mesh, triangulate, cellIdLabels);
}


Foam::cuttingPlane::cuttingPlane
(
    const plane& pln,
    const primitiveMesh& mesh,
    const bool triangulate,
    const labelUList& cellIdLabels
)
:
    plane(pln)
{
    performCut(mesh, triangulate, cellIdLabels);
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::cuttingPlane::performCut
(
    const primitiveMesh& mesh,
    const bool triangulate,
    bitSet&& cellIdLabels
)
{
    MeshStorage::clear();
    meshCells_.clear();

    // Pre-populate with restriction
    bitSet cellCuts(std::move(cellIdLabels));

    if (cellCuts.size())
    {
        cellCuts.resize(mesh.nCells());
    }

    // Classification of points with respect to the plane
    PackedList<2> sides = classifySides(*this, mesh.points());

    // Determine cells that are (likely) cut
    // - some ambiguity when plane is exactly between cells
    const label nFaceCuts = calcCellCuts(mesh, sides, cellCuts);

    // Find closed loop from cell cuts
    walkCellCuts(mesh, cellCuts, sides, triangulate, nFaceCuts);
}


void Foam::cuttingPlane::performCut
(
    const primitiveMesh& mesh,
    const bool triangulate,
    const bitSet& cellIdLabels
)
{
    bitSet subsetCells(cellIdLabels);

    performCut(mesh, triangulate, std::move(subsetCells));
}


void Foam::cuttingPlane::performCut
(
    const primitiveMesh& mesh,
    const bool triangulate,
    const labelUList& cellIdLabels
)
{
    bitSet subsetCells;

    if (notNull(cellIdLabels))
    {
        // Pre-populate with restriction
        subsetCells.resize(mesh.nCells());
        subsetCells.set(cellIdLabels);
    }

    performCut(mesh, triangulate, std::move(subsetCells));
}


void Foam::cuttingPlane::remapFaces(const labelUList& faceMap)
{
    if (notNull(faceMap) && !faceMap.empty())
    {
        MeshStorage::remapFaces(faceMap);

        List<label> newCutCells(faceMap.size());
        forAll(faceMap, facei)
        {
            newCutCells[facei] = meshCells_[faceMap[facei]];
        }
        meshCells_.transfer(newCutCells);
    }
}


// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //

void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
{
    if (this == &rhs)
    {
        FatalErrorInFunction
            << "Attempted assignment to self"
            << abort(FatalError);
    }

    static_cast<MeshStorage&>(*this) = rhs;
    static_cast<plane&>(*this) = rhs;
    meshCells_ = rhs.meshCells();
}


// ************************************************************************* //