polyMeshZipUpCells.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-2016 OpenFOAM Foundation
    Copyright (C) 2022 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 "polyMeshZipUpCells.H"
#include "polyMesh.H"
#include "Time.H"
#include "CircularBuffer.H"
#include "DynamicList.H"

// #define DEBUG_ZIPUP 1
// #define DEBUG_CHAIN 1
// #define DEBUG_ORDER 1

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

bool Foam::polyMeshZipUpCells(polyMesh& mesh)
{
    if (polyMesh::debug)
    {
        Info<< "bool polyMeshZipUpCells(polyMesh& mesh) const: "
            << "zipping up topologically open cells" << endl;
    }

    // Algorithm:
    // Take the original mesh and visit all cells.  For every cell
    // calculate the edges of all faces on the cells.  A cell is
    // correctly topologically closed when all the edges are referenced
    // by exactly two faces.  If the edges are referenced only by a
    // single face, additional vertices need to be inserted into some
    // of the faces (topological closedness).  If an edge is
    // referenced by more that two faces, there is an error in
    // topological closedness.
    // Point insertion into the faces is done by attempting to create
    // closed loops and inserting the intermediate points into the
    // defining edge
    // Note:
    // The algorithm is recursive and changes the mesh faces in each
    // pass.  It is therefore essential to discard the addressing
    // after every pass.  The algorithm is completed when the mesh
    // stops changing.

    label nChangedFacesInMesh = 0;
    label nCycles = 0;

    labelHashSet problemCells;

    DynamicList<bool> singleEdgeUsage(256);
    CircularBuffer<label> pointChain(256);

    do
    {
        nChangedFacesInMesh = 0;

        const cellList& Cells = mesh.cells();
        const pointField& Points = mesh.points();

        faceList newFaces = mesh.faces();

        const faceList& oldFaces = mesh.faces();
        const labelListList& pFaces = mesh.pointFaces();

        forAll(Cells, celli)
        {
            const labelList& curFaces = Cells[celli];
            const edgeList cellEdges = Cells[celli].edges(oldFaces);
            const labelList cellPoints = Cells[celli].labels(oldFaces);

            // Find the edges used only once in the cell

            labelList edgeUsage(cellEdges.size(), Zero);

            forAll(curFaces, facei)
            {
                edgeList curFaceEdges = oldFaces[curFaces[facei]].edges();

                forAll(curFaceEdges, faceEdgeI)
                {
                    const edge& curEdge = curFaceEdges[faceEdgeI];

                    forAll(cellEdges, cellEdgeI)
                    {
                        if (cellEdges[cellEdgeI] == curEdge)
                        {
                            edgeUsage[cellEdgeI]++;
                            break;
                        }
                    }
                }
            }

            edgeList singleEdges(cellEdges.size());
            label nSingleEdges = 0;

            forAll(edgeUsage, edgeI)
            {
                if (edgeUsage[edgeI] == 1)
                {
                    singleEdges[nSingleEdges] = cellEdges[edgeI];
                    nSingleEdges++;
                }
                else if (edgeUsage[edgeI] != 2)
                {
                    WarningInFunction
                        << "edge " << cellEdges[edgeI] << " in cell " << celli
                        << " used " << edgeUsage[edgeI] << " times. " << nl
                        << "Should be 1 or 2 - serious error "
                        << "in mesh structure. " << endl;

                    #ifdef DEBUG_ZIPUP
                    forAll(curFaces, facei)
                    {
                        Info<< "face: " << oldFaces[curFaces[facei]]
                            << endl;
                    }

                    Info<< "Cell edges: " << cellEdges << nl
                        << "Edge usage: " << edgeUsage << nl
                        << "Cell points: " << cellPoints << endl;

                    forAll(cellPoints, cpI)
                    {
                        Info<< "vertex create \"" << cellPoints[cpI]
                            << "\" coordinates "
                            << Points[cellPoints[cpI]] << endl;
                    }
                    #endif

                    // Gather the problem cell
                    problemCells.insert(celli);
                }
            }

            // Check if the cell is already zipped up
            if (nSingleEdges == 0) continue;

            singleEdges.setSize(nSingleEdges);

            #ifdef DEBUG_ZIPUP
            Info<< "Cell " << celli << endl;

            forAll(curFaces, facei)
            {
                Info<< "face: " << oldFaces[curFaces[facei]] << endl;
            }

            Info<< "Cell edges: " << cellEdges << nl
                << "Edge usage: " << edgeUsage << nl
                << "Single edges: " << singleEdges << nl
                << "Cell points: " << cellPoints << endl;

            forAll(cellPoints, cpI)
            {
                Info<< "vertex create \"" << cellPoints[cpI]
                    << "\" coordinates "
                    << points()[cellPoints[cpI]] << endl;
            }
            #endif

            // Loop through all single edges and mark the points they use
            // points marked twice are internal to edge; those marked more than
            // twice are corners

            labelList pointUsage(cellPoints.size(), Zero);

            forAll(singleEdges, edgeI)
            {
                const edge& curEdge = singleEdges[edgeI];

                forAll(cellPoints, pointi)
                {
                    if
                    (
                        cellPoints[pointi] == curEdge.start()
                     || cellPoints[pointi] == curEdge.end()
                    )
                    {
                        pointUsage[pointi]++;
                    }
                }
            }

            singleEdgeUsage.resize_nocopy(singleEdges.size());
            singleEdgeUsage = false;

            // loop through all edges and eliminate the ones that are
            // blocked out
            forAll(singleEdges, edgeI)
            {
                bool blockedHead = false;
                bool blockedTail = false;

                label newEdgeStart = singleEdges[edgeI].start();
                label newEdgeEnd = singleEdges[edgeI].end();

                // check that the edge has not got all ends blocked
                forAll(cellPoints, pointi)
                {
                    if (cellPoints[pointi] == newEdgeStart)
                    {
                        if (pointUsage[pointi] > 2)
                        {
                            blockedHead = true;
                        }
                    }
                    else if (cellPoints[pointi] == newEdgeEnd)
                    {
                        if (pointUsage[pointi] > 2)
                        {
                            blockedTail = true;
                        }
                    }
                }

                if (blockedHead && blockedTail)
                {
                    // Eliminating edge singleEdges[edgeI] as blocked
                    singleEdgeUsage[edgeI] = true;
                }
            }

            // Go through the points and start from the point used twice
            // check all the edges to find the edges starting from this point
            // add the

            labelListList edgesToInsert(singleEdges.size());
            label nEdgesToInsert = 0;

            // Find a good edge
            forAll(singleEdges, edgeI)
            {
                pointChain.clear();

                bool blockHead = false;
                bool blockTail = false;

                if (!singleEdgeUsage[edgeI])
                {
                    // found a new edge
                    singleEdgeUsage[edgeI] = true;

                    label newEdgeStart = singleEdges[edgeI].start();
                    label newEdgeEnd = singleEdges[edgeI].end();

                    pointChain.push_front(newEdgeStart);
                    pointChain.push_back(newEdgeEnd);

                    #ifdef DEBUG_CHAIN
                    Info<< "found edge to start with: "
                        << singleEdges[edgeI] << endl;
                    #endif

                    // Check if head or tail are blocked
                    forAll(cellPoints, pointi)
                    {
                        if (cellPoints[pointi] == newEdgeStart)
                        {
                            if (pointUsage[pointi] > 2)
                            {
                                #ifdef DEBUG_CHAIN
                                Info<< "start head blocked" << endl;
                                #endif

                                blockHead = true;
                            }
                        }
                        else if (cellPoints[pointi] == newEdgeEnd)
                        {
                            if (pointUsage[pointi] > 2)
                            {
                                #ifdef DEBUG_CHAIN
                                Info<< "start tail blocked" << endl;
                                #endif

                                blockTail = true;
                            }
                        }
                    }

                    bool stopSearching = false;

                    // Go through the unused edges and try to chain them up
                    do
                    {
                        stopSearching = false;

                        forAll(singleEdges, addEdgeI)
                        {
                            if (!singleEdgeUsage[addEdgeI])
                            {
                                // Grab start and end of the candidate
                                label addStart =
                                    singleEdges[addEdgeI].start();

                                label addEnd =
                                    singleEdges[addEdgeI].end();

                                #ifdef DEBUG_CHAIN
                                Info<< "Trying candidate "
                                    << singleEdges[addEdgeI] << endl;
                                #endif

                                // Try to add the edge onto the head
                                if (!blockHead)
                                {
                                    if (pointChain.front() == addStart)
                                    {
                                        // Added at start mark as used
                                        pointChain.push_front(addEnd);

                                        singleEdgeUsage[addEdgeI] = true;
                                    }
                                    else if (pointChain.front() == addEnd)
                                    {
                                        pointChain.push_front(addStart);

                                        singleEdgeUsage[addEdgeI] = true;
                                    }
                                }

                                // Try the other end only if the first end
                                // did not add it
                                if (!blockTail && !singleEdgeUsage[addEdgeI])
                                {
                                    if (pointChain.back() == addStart)
                                    {
                                        // Added at start mark as used
                                        pointChain.push_back(addEnd);

                                        singleEdgeUsage[addEdgeI] = true;
                                    }
                                    else if (pointChain.back() == addEnd)
                                    {
                                        pointChain.push_back(addStart);

                                        singleEdgeUsage[addEdgeI] = true;
                                    }
                                }

                                // check if the new head or tail are blocked
                                label curEdgeStart = pointChain.front();
                                label curEdgeEnd = pointChain.back();

                                #ifdef DEBUG_CHAIN
                                Info<< "curEdgeStart: " << curEdgeStart
                                    << " curEdgeEnd: " << curEdgeEnd << endl;
                                #endif

                                forAll(cellPoints, pointi)
                                {
                                    if (cellPoints[pointi] == curEdgeStart)
                                    {
                                        if (pointUsage[pointi] > 2)
                                        {
                                            #ifdef DEBUG_CHAIN
                                            Info<< "head blocked" << endl;
                                            #endif

                                            blockHead = true;
                                        }
                                    }
                                    else if (cellPoints[pointi] == curEdgeEnd)
                                    {
                                        if (pointUsage[pointi] > 2)
                                        {
                                            #ifdef DEBUG_CHAIN
                                            Info<< "tail blocked" << endl;
                                            #endif

                                            blockTail = true;
                                        }
                                    }
                                }

                                // Check if the loop is closed
                                if (curEdgeStart == curEdgeEnd)
                                {
                                    #ifdef DEBUG_CHAIN
                                    Info<< "closed loop" << endl;
                                    #endif

                                    pointChain.pop_front();

                                    blockHead = true;
                                    blockTail = true;

                                    stopSearching = true;
                                }

                                #ifdef DEBUG_CHAIN
                                Info<< "current pointChain: "
                                    << pointChain << endl;
                                #endif

                                if (stopSearching) break;
                            }
                        }
                    } while (stopSearching);
                }

                #ifdef DEBUG_CHAIN
                Info<< "completed patch chain: " << pointChain << endl;
                #endif

                if (pointChain.size() > 2)
                {
                    edgesToInsert[nEdgesToInsert] = pointChain.list();
                    nEdgesToInsert++;
                }
            }

            edgesToInsert.setSize(nEdgesToInsert);

            #ifdef DEBUG_ZIPUP
            Info<< "edgesToInsert: " << edgesToInsert << endl;
            #endif

            // Insert the edges into a list of faces
            forAll(edgesToInsert, edgeToInsertI)
            {
                // Order the points of the edge
                // Warning: the ordering must be parametric, because in
                // the case of multiple point insertion onto the same edge
                // it is possible to get non-cyclic loops
                //

                const labelList& unorderedEdge = edgesToInsert[edgeToInsertI];

                scalarField dist(unorderedEdge.size());

                // Calculate distance
                point startPoint = Points[unorderedEdge[0]];
                dist[0] = 0;

                vector dir = Points[unorderedEdge.last()] - startPoint;

                for (label i = 1; i < dist.size(); i++)
                {
                    dist[i] = (Points[unorderedEdge[i]] - startPoint) & dir;
                }

                // Sort points
                labelList orderedEdge(unorderedEdge.size(), -1);
                boolList used(unorderedEdge.size(), false);

                forAll(orderedEdge, epI)
                {
                    label nextPoint = -1;
                    scalar minDist = GREAT;

                    forAll(dist, i)
                    {
                        if (!used[i] && dist[i] < minDist)
                        {
                            minDist = dist[i];
                            nextPoint = i;
                        }
                    }

                    // Insert the next point
                    orderedEdge[epI] = unorderedEdge[nextPoint];
                    used[nextPoint] = true;
                }

                #ifdef DEBUG_ORDER
                Info<< "unorderedEdge: " << unorderedEdge << nl
                    << "orderedEdge: " << orderedEdge << endl;
                #endif

                // check for duplicate points in the ordered edge
                forAll(orderedEdge, checkI)
                {
                    for
                    (
                        label checkJ = checkI + 1;
                        checkJ < orderedEdge.size();
                        checkJ++
                    )
                    {
                        if (orderedEdge[checkI] == orderedEdge[checkJ])
                        {
                            WarningInFunction
                                << "Duplicate point found in edge to insert. "
                                << nl << "Point: " << orderedEdge[checkI]
                                << " edge: " << orderedEdge << endl;

                            problemCells.insert(celli);
                        }
                    }
                }

                edge testEdge
                (
                    orderedEdge[0],
                    orderedEdge.last()
                );

                // In order to avoid edge-to-edge comparison, get faces using
                // point-face addressing in two goes.
                const labelList& startPF = pFaces[testEdge.start()];
                const labelList& endPF = pFaces[testEdge.start()];

                labelList facesSharingEdge(startPF.size() + endPF.size());
                label nfse = 0;

                forAll(startPF, pfI)
                {
                    facesSharingEdge[nfse++] = startPF[pfI];
                }
                forAll(endPF, pfI)
                {
                    facesSharingEdge[nfse++] = endPF[pfI];
                }

                forAll(facesSharingEdge, facei)
                {
                    bool faceChanges = false;

                    // Label of the face being analysed
                    const label currentFaceIndex = facesSharingEdge[facei];

                    const edgeList curFaceEdges =
                        oldFaces[currentFaceIndex].edges();

                    forAll(curFaceEdges, cfeI)
                    {
                        if (curFaceEdges[cfeI] == testEdge)
                        {
                            faceChanges = true;
                            break;
                        }
                    }

                    if (faceChanges)
                    {
                        nChangedFacesInMesh++;
                        // In order to avoid losing point from multiple
                        // insertions into the same face, the new face
                        // will be change incrementally.
                        // 1) Check if all the internal points of the edge
                        // to add already exist in the face. If so, the
                        // edge has already been included 2) Check if the
                        // point insertion occurs on an edge which is
                        // still untouched. If so, simply insert
                        // additional points into the face.  3) If not,
                        // the edge insertion occurs on an already
                        // modified edge. ???

                        face& newFace = newFaces[currentFaceIndex];

                        bool allPointsPresent = true;

                        forAll(orderedEdge, oeI)
                        {
                            bool curPointFound = false;

                            forAll(newFace, nfI)
                            {
                                if (newFace[nfI] == orderedEdge[oeI])
                                {
                                    curPointFound = true;
                                    break;
                                }
                            }

                            allPointsPresent =
                                allPointsPresent && curPointFound;
                        }

                        #ifdef DEBUG_ZIPUP
                        if (allPointsPresent)
                        {
                            Info<< "All points present" << endl;
                        }
                        #endif

                        if (!allPointsPresent)
                        {
                            // Not all points are already present.  The
                            // new edge will need to be inserted into the
                            // face.

                            // Check to see if a new edge fits onto an
                            // untouched edge of the face.  Make sure the
                            // edges are grabbed before the face is
                            // resized.
                            edgeList newFaceEdges = newFace.edges();

                            #ifdef DEBUG_ZIPUP
                            Info<< "Not all points present." << endl;
                            #endif

                            label nNewFacePoints = 0;

                            bool edgeAdded = false;

                            forAll(newFaceEdges, curFacEdgI)
                            {
                                // Does the current edge change?
                                if (newFaceEdges[curFacEdgI] == testEdge)
                                {
                                    // Found an edge match
                                    edgeAdded = true;

                                    // Resize the face to accept additional
                                    // points
                                    newFace.setSize
                                    (
                                        newFace.size()
                                      + orderedEdge.size() - 2
                                    );

                                    if
                                    (
                                        newFaceEdges[curFacEdgI].start()
                                     == testEdge.start()
                                    )
                                    {
                                        // insertion in ascending order
                                        for
                                        (
                                            label i = 0;
                                            i < orderedEdge.size() - 1;
                                            i++
                                        )
                                        {
                                            newFace[nNewFacePoints] =
                                                orderedEdge[i];
                                            nNewFacePoints++;
                                        }
                                    }
                                    else
                                    {
                                        // insertion in reverse order
                                        for
                                        (
                                            label i = orderedEdge.size() - 1;
                                            i > 0;
                                            i--
                                        )
                                        {
                                            newFace[nNewFacePoints] =
                                                orderedEdge[i];
                                            nNewFacePoints++;
                                        }
                                    }
                                }
                                else
                                {
                                    // Does not fit onto this edge.
                                    // Copy the next point into the face
                                    newFace[nNewFacePoints] =
                                        newFaceEdges[curFacEdgI].start();
                                    nNewFacePoints++;
                                }
                            }

                            #ifdef DEBUG_ZIPUP
                            Info<< "oldFace: "
                                << oldFaces[currentFaceIndex] << nl
                                << "newFace: " << newFace << endl;
                            #endif

                            // Check for duplicate points in the new face
                            forAll(newFace, checkI)
                            {
                                for
                                (
                                    label checkJ = checkI + 1;
                                    checkJ < newFace.size();
                                    checkJ++
                                )
                                {
                                    if (newFace[checkI] == newFace[checkJ])
                                    {
                                        WarningInFunction
                                            << "Duplicate point found "
                                            << "in the new face. " << nl
                                            << "Point: "
                                            << orderedEdge[checkI]
                                            << " face: "
                                            << newFace << endl;

                                        problemCells.insert(celli);
                                    }
                                }
                            }

                            // Check if the edge is added.
                            // If not, then it comes on top of an already
                            // modified edge and they need to be
                            // merged in together.
                            if (!edgeAdded)
                            {
                                Info<< "This edge modifies an already modified "
                                    << "edge.  Point insertions skipped."
                                    << endl;
                            }
                        }
                    }
                }
            }
        }

        if (problemCells.size())
        {
            // This cycle has failed.  Print out the problem cells
            FatalErrorInFunction
                << "Found " << problemCells.size() << " problem cells." << nl
                << "Cells: " << problemCells.sortedToc()
                << abort(FatalError);
        }

        Info<< "Cycle " << ++nCycles
            << " changed " << nChangedFacesInMesh << " faces." << endl;


        const polyBoundaryMesh& bMesh = mesh.boundaryMesh();

        // Reset the polyMesh. Number of points/faces/cells/patches stays the
        // same, only the faces themselves have changed so clear all derived
        // (edge, point) addressing.

        // Collect the patch sizes
        labelList patchSizes(bMesh.size(), Zero);
        labelList patchStarts(bMesh.size(), Zero);

        forAll(bMesh, patchi)
        {
            patchSizes[patchi] = bMesh[patchi].size();
            patchStarts[patchi] = bMesh[patchi].start();
        }

        // Reset the mesh. Number of active faces is one beyond the last patch
        // (patches guaranteed to be in increasing order)
        mesh.resetPrimitives
        (
            autoPtr<pointField>(),  // <- null: leaves points untouched
            autoPtr<faceList>::New(std::move(newFaces)),
            autoPtr<labelList>(),   // <- null: leaves owner untouched
            autoPtr<labelList>(),   // <- null: leaves neighbour untouched
            patchSizes,
            patchStarts,
            true                // boundary forms valid boundary mesh.
        );

        // Reset any addressing on face zones.
        mesh.faceZones().clearAddressing();

        // Clear the addressing
        mesh.clearOut();

    } while (nChangedFacesInMesh > 0 || nCycles > 100);

    // Flags the mesh files as being changed
    mesh.setInstance(mesh.time().timeName());

    if (nChangedFacesInMesh > 0)
    {
        FatalErrorInFunction
            << "with the original mesh"
            << abort(FatalError);
    }

    return nCycles != 1;
}


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