removeCellLayer.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) 2018-2024 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 "layerAdditionRemoval.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "polyTopoChange.H"
#include "oppositeFace.H"
#include "polyTopoChanger.H"
#include "polyRemoveCell.H"
#include "polyRemoveFace.H"
#include "polyRemovePoint.H"
#include "polyModifyFace.H"
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
bool Foam::layerAdditionRemoval::validCollapse() const
{
    // Check for valid layer collapse
    // - no boundary-to-boundary collapse
 
    if (debug)
    {
        Pout<< "Checking layer collapse for object " << name() << endl;
    }
 
    // Grab the face collapse mapping
    const polyMesh& mesh = topoChanger().mesh();
 
    const labelList& ftc = facesPairing();
    const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
 
    label nBoundaryHits = 0;
 
    forAll(mf, facei)
    {
        if
        (
            !mesh.isInternalFace(mf[facei])
         && !mesh.isInternalFace(ftc[facei])
        )
        {
            nBoundaryHits++;
        }
    }
 
 
    if (debug)
    {
        Pout<< "Finished checking layer collapse for object "
            << name() <<".  Number of boundary-on-boundary hits: "
            << nBoundaryHits << endl;
    }
 
    return !returnReduceOr(nBoundaryHits);
}
 
 
void Foam::layerAdditionRemoval::removeCellLayer
(
    polyTopoChange& ref
) const
{
    // Algorithm for layer removal.  Second phase: topological change
    // 2)  Go through all the faces of the master cell layer and remove
    //     the ones that are not in the master face zone.
    //
    // 3)  Grab all the faces coming out of points that are collapsed
    //     and select the ones that are not marked for removal.  Go
    //     through the remaining faces and replace the point to remove by
    //     the equivalent point in the master face zone.
    if (debug)
    {
        Pout<< "Removing the cell layer for object " << name() << endl;
    }
 
    const polyMesh& mesh = topoChanger().mesh();
 
    const labelList& own = mesh.faceOwner();
    const labelList& nei = mesh.faceNeighbour();
 
    const labelList& ptc = pointsPairing();
    const labelList& ftc = facesPairing();
 
    // Remove all the cells from the master layer
    const labelList& mc =
        mesh.faceZones()[faceZoneID_.index()].masterCells();
 
    forAll(mc, facei)
    {
        label slaveSideCell = own[ftc[facei]];
 
        if (mesh.isInternalFace(ftc[facei]) && slaveSideCell == mc[facei])
        {
            // Owner cell of the face is being removed.
            // Grab the neighbour instead
            slaveSideCell = nei[ftc[facei]];
        }
 
        ref.setAction(polyRemoveCell(mc[facei], slaveSideCell));
    }
 
    // Remove all the faces from the master layer cells which are not in
    // the master face layer
    labelHashSet facesToRemoveMap(mc.size()*primitiveMesh::facesPerCell_);
 
    const cellList& cells = mesh.cells();
 
    forAll(mc, celli)
    {
        const cell& curCell = cells[mc[celli]];
 
        forAll(curCell, facei)
        {
            // Check if the face is in the master zone.  If not, remove it
            if
            (
                mesh.faceZones().whichZone(curCell[facei])
             != faceZoneID_.index()
            )
            {
                facesToRemoveMap.insert(curCell[facei]);
            }
        }
    }
 
    for (const label facei : facesToRemoveMap)
    {
        ref.setAction(polyRemoveFace(facei));
    }
 
    // Remove all points that will be collapsed
    for (const label pointi : ptc)
    {
        ref.setAction(polyRemovePoint(pointi));
    }
 
    // Grab all faces coming off points to be deleted.  If the face
    // has not been deleted, replace the removed point with the
    // equivalent from the master layer.
 
    // Make a map of all point to be removed, giving the master point label
    // for each of them
 
    const labelList& meshPoints =
        mesh.faceZones()[faceZoneID_.index()].patch().meshPoints();
 
    Map<label> removedPointMap(ptc, meshPoints);
 
    const labelListList& pf = mesh.pointFaces();
 
    const faceList& faces = mesh.faces();
 
    // Make a list of faces to be modified using the map to avoid duplicates
    labelHashSet facesToModify(ptc.size()*primitiveMesh::facesPerPoint_);
 
    forAll(ptc, pointi)
    {
        const labelList& curFaces = pf[ptc[pointi]];
 
        forAll(curFaces, facei)
        {
            if (!facesToRemoveMap.found(curFaces[facei]))
            {
                facesToModify.insert(curFaces[facei]);
            }
        }
    }
 
    labelList ftm = facesToModify.toc();
 
    if (debug > 1)
    {
        Pout<< "faces to modify: " << ftm << endl;
    }
 
    forAll(ftm, facei)
    {
        // For every face to modify, copy the face and re-map the vertices.
        // It is known all the faces will be changed since they hang off
        // re-mapped vertices
        label curFaceID = ftm[facei];
 
        face newFace(faces[curFaceID]);
 
        forAll(newFace, pointi)
        {
            const auto rpmIter = removedPointMap.cfind(newFace[pointi]);
 
            if (rpmIter.good())
            {
                // Point mapped. Replace it
                newFace[pointi] = rpmIter();
            }
        }
 
        if (debug > 1)
        {
            Pout<< "face label: " << curFaceID
                << " old face: " << faces[curFaceID]
                << " new face: " << newFace << endl;
        }
 
        // Get face zone and its flip
        label modifiedFaceZone = mesh.faceZones().whichZone(curFaceID);
        bool modifiedFaceZoneFlip = false;
 
        if (modifiedFaceZone >= 0)
        {
            const faceZone& fz = mesh.faceZones()[modifiedFaceZone];
            modifiedFaceZoneFlip = fz.flipMap()[fz.whichFace(curFaceID)];
        }
 
        label newNeighbour = -1;
 
        if (curFaceID < mesh.nInternalFaces())
        {
            newNeighbour = nei[curFaceID];
        }
 
        // Modify the face
        ref.setAction
        (
            polyModifyFace
            (
                newFace,                // modified face
                curFaceID,              // label of face being modified
                own[curFaceID],         // owner
                newNeighbour,           // neighbour
                false,                  // face flip
                mesh.boundaryMesh().whichPatch(curFaceID),// patch for face
                false,                  // remove from zone
                modifiedFaceZone,       // zone for face
                modifiedFaceZoneFlip    // face flip in zone
            )
        );
    }
 
    // Modify the faces in the master layer to point past the removed cells
 
    const labelList& mf = mesh.faceZones()[faceZoneID_.index()];
    const boolList& mfFlip = mesh.faceZones()[faceZoneID_.index()].flipMap();
 
    forAll(mf, facei)
    {
        // Grab the owner and neighbour of the faces to be collapsed and get rid
        // of the cell to be removed
        label masterSideCell = own[mf[facei]];
 
        if (masterSideCell == mc[facei])
        {
            if (mesh.isInternalFace(mf[facei]))
            {
                // Owner cell of the face is being removed.
                // Grab the neighbour instead
                masterSideCell = nei[mf[facei]];
            }
            else
            {
                masterSideCell = -1;
            }
        }
 
        label slaveSideCell = own[ftc[facei]];
 
        if (slaveSideCell == mc[facei])
        {
            if (mesh.isInternalFace(ftc[facei]))
            {
                // Owner cell of the face is being removed.
                // Grab the neighbour instead
                slaveSideCell = nei[ftc[facei]];
            }
            else
            {
                slaveSideCell = -1;
            }
        }
 
        // Find out if the face needs to be flipped
        label newOwner = -1;
        label newNeighbour = -1;
        bool flipFace = false;
        label newPatchID = -1;
        label newZoneID = -1;
 
        // Is any of the faces a boundary face?  If so, grab the patch
        // A boundary-to-boundary collapse is checked for in validCollapse()
        // and cannot happen here.
 
        if (!mesh.isInternalFace(mf[facei]))
        {
            // Master is the boundary face: it gets a new owner but no flip
            newOwner = slaveSideCell;
            newNeighbour = -1;
            flipFace = false;
            newPatchID = mesh.boundaryMesh().whichPatch(mf[facei]);
            newZoneID = mesh.faceZones().whichZone(mf[facei]);
        }
        else if (!mesh.isInternalFace(ftc[facei]))
        {
            // Slave is the boundary face: grab its patch
            newOwner = slaveSideCell;
            newNeighbour = -1;
 
            // Find out if the face flip is necessary
            if (own[mf[facei]] == slaveSideCell)
            {
                flipFace = false;
            }
            else
            {
                flipFace = true;
            }
 
            newPatchID = mesh.boundaryMesh().whichPatch(ftc[facei]);
 
            // The zone of the master face is preserved
            newZoneID = mesh.faceZones().whichZone(mf[facei]);
        }
        else
        {
            // Both faces are internal: flip is decided based on which of the
            // new cells around it has a lower label
            newOwner = min(masterSideCell, slaveSideCell);
            newNeighbour = max(masterSideCell, slaveSideCell);
 
            if (newOwner == own[mf[facei]] || newNeighbour == nei[mf[facei]])
            {
                flipFace = false;
            }
            else
            {
                flipFace = true;
            }
 
            newPatchID = -1;
 
            // The zone of the master face is preserved
            newZoneID = mesh.faceZones().whichZone(mf[facei]);
        }
 
        // Modify the face and flip if necessary
        face newFace = faces[mf[facei]];
        bool zoneFlip = mfFlip[facei];
 
        if (flipFace)
        {
            newFace.flip();
            zoneFlip = !zoneFlip;
        }
 
        if (debug > 1)
        {
            Pout<< "Modifying face " << mf[facei]
                << " newFace: " << newFace << nl
                << " newOwner: " << newOwner
                << " newNeighbour: " << newNeighbour
                << " flipFace: " << flipFace
                << " newPatchID: " << newPatchID
                << " newZoneID: " << newZoneID << nl
                << " oldOwn: " << own[mf[facei]];
            if (newPatchID == -1)
            {
                Pout<< " oldNei: " << nei[mf[facei]];
            }
            Pout<< endl;
        }
 
        ref.setAction
        (
            polyModifyFace
            (
                newFace,       // modified face
                mf[facei],     // label of face being modified
                newOwner,      // owner
                newNeighbour,  // neighbour
                flipFace,      // flip
                newPatchID,    // patch for face
                false,         // remove from zone
                newZoneID,     // new zone
                zoneFlip       // face flip in zone
            )
        );
    }
}
 
 
// ************************************************************************* //