createPolyCells.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) 2016 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/>.
 
Description
    create cellPolys
    - use pointCells when searching for connectivity
    - initialize the cell connectivity with '-1'
    - find both cell faces corresponding to the baffles and mark them
      to prevent a connection
    - standard connectivity checks
 
    - added baffle support
 
\*---------------------------------------------------------------------------*/
 
#include "meshReader.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
void Foam::meshReader::createPolyCells()
{
    // loop through all cell faces and create connectivity. This will produce
    // a global face list and will describe all cells as lists of face labels
 
    const faceListList& cFaces = cellFaces();
 
    // count the maximum number of faces and set the size of the cellPolys_
    cellPolys_.setSize(cFaces.size());
 
    label maxFaces = 0;
 
    forAll(cellPolys_, celli)
    {
        cellPolys_[celli].setSize(cFaces[celli].size(), -1);
 
        maxFaces += cFaces[celli].size();
    }
 
    Info<< "Maximum possible number of faces in mesh: " << maxFaces << endl;
 
    meshFaces_.setSize(maxFaces);
 
    // set reference to point-cell addressing
    const labelListList& ptCells = pointCells();
 
    // size the baffle lists and initialize to -1
    baffleIds_.setSize(baffleFaces_.size());
    forAll(baffleIds_, baffleI)
    {
        baffleIds_[baffleI].setSize(2);
    }
 
    // block off baffles first
    //
    // To prevent internal faces, we'll mark the cell faces
    // with negative cell ids (offset by nCells).
    // eg,
    //    celli = -(nCells + baffleI)
    //
    // To distinguish these from the normal '-1' marker, we require
    //    celli = -(nCells + baffleI) < -1
    //
    // This condition is met provided that nCells > 1.
    // ie., baffles require at least 2 volume cells
 
    label baffleOffset = cFaces.size();
    forAll(baffleFaces_, baffleI)
    {
        label celli = -(baffleOffset + baffleI);
        const face& curFace = baffleFaces_[baffleI];
 
        // get the list of labels
        const labelList& curPoints = curFace;
 
        // a baffle is a single face - only need to match one face
        // get the list of cells sharing this point
        const labelList& curNeighbours = ptCells[curPoints[0]];
 
        label nNeighbours = 0;
 
        // For all neighbours
        forAll(curNeighbours, neiI)
        {
            label curNei = curNeighbours[neiI];
 
            // get the list of search faces
            const faceList& searchFaces = cFaces[curNei];
 
            forAll(searchFaces, neiFacei)
            {
                int cmp = face::compare(curFace, searchFaces[neiFacei]);
 
                if (cmp)
                {
                    // maintain baffle orientation
                    // side0: baffle normal same as attached face
                    // side1: baffle normal opposite from attached face
                    //
                    label side = 0;
                    if (cmp < 0)
                    {
                        side = 1;
                    }
 
#ifdef DEBUG_FACE_ORDERING
                    Info<< "cmp " << cmp << " matched " << curFace
                        << " with " << searchFaces[neiFacei]
                        << endl;
 
 
                    Info<< "match " << baffleI
                        << " (" << origCellId_[baffleOffset+baffleI] << ")"
                        << " side " << side
                        << " against cell " << curNei
                        << " face " << neiFacei
                        << " curFace " << curFace[1]
                        << " neiFace " << searchFaces[neiFacei][1]
                        << endl;
#endif
 
                    if (baffleIds_[baffleI][side].notUsed())
                    {
                        baffleIds_[baffleI][side] = cellFaceIdentifier
                        (
                            curNei,
                            neiFacei
                        );
 
                        nNeighbours++;
                    }
                    else
                    {
                        Info<< "multiple matches for side " << side
                            << " of baffle " << baffleI
                            << " (original cell "
                            << origCellId_[baffleOffset+baffleI] << ")"
                            << endl;
                    }
                    break;
                }
            }
            if (nNeighbours >= 2) break;
        }
 
        if (nNeighbours == 2)
        {
            for (label side = 0; side < nNeighbours; ++side)
            {
                label neiCell = baffleIds_[baffleI][side].cellId();
                label neiFace = baffleIds_[baffleI][side].faceId();
 
                if (baffleIds_[baffleI][side].used())
                {
                    cellPolys_[neiCell][neiFace] = celli;
                }
            }
        }
        else
        {
            Info<< "drop baffle " << baffleI
                << " (original cell "
                << origCellId_[baffleOffset+baffleI] << ")"
                << " with " << nNeighbours << " neighbours" << endl;
 
            baffleFaces_[baffleI].clear();
            baffleIds_[baffleI].clear();
        }
    }
 
#ifdef DEBUG_CELLPOLY
    Info<< "cellPolys_" << cellPolys_ << endl;
    Info<< "baffleFaces_" << baffleFaces_ << endl;
    Info<< "baffleIds_"   << baffleIds_ << endl;
#endif
 
    bool found = false;
 
    nInternalFaces_ = 0;
 
    forAll(cFaces, celli)
    {
        // Note:
        // Insertion cannot be done in one go as the faces need to be
        // added into the list in the increasing order of neighbour
        // cells.  Therefore, all neighbours will be detected first
        // and then added in the correct order.
 
        const faceList& curFaces = cFaces[celli];
 
        // Record the neighbour cell
        labelList neiCells(curFaces.size(), -1);
 
        // Record the face of neighbour cell
        labelList faceOfNeiCell(curFaces.size(), -1);
 
        label nNeighbours = 0;
 
        // For all faces ...
        forAll(curFaces, facei)
        {
            // Skip already matched faces or those tagged by baffles
            if (cellPolys_[celli][facei] != -1) continue;
 
            found = false;
 
            const face& curFace = curFaces[facei];
 
            // get the list of labels
            const labelList& curPoints = curFace;
 
            // For all points
            forAll(curPoints, pointi)
            {
                // get the list of cells sharing this point
                const labelList& curNeighbours = ptCells[curPoints[pointi]];
 
                // For all neighbours
                forAll(curNeighbours, neiI)
                {
                    label curNei = curNeighbours[neiI];
 
                    // reject neighbours with the lower label. This should
                    // also reject current cell.
                    if (curNei > celli)
                    {
                        // get the list of search faces
                        const faceList& searchFaces = cFaces[curNei];
 
                        forAll(searchFaces, neiFacei)
                        {
                            if (searchFaces[neiFacei] == curFace)
                            {
                                // Record the neighbour cell and face
                                neiCells[facei] = curNei;
                                faceOfNeiCell[facei] = neiFacei;
                                nNeighbours++;
#ifdef DEBUG_FACE_ORDERING
                                Info<< " cell " << celli
                                    << " face " << facei
                                    << " point " << pointi
                                    << " nei " << curNei
                                    << " neiFace " << neiFacei
                                    << endl;
#endif
                                found = true;
                                break;
                            }
                        }
                        if (found) break;
                    }
                    if (found) break;
                }
                if (found) break;
            } // End of current points
        } // End of current faces
 
        // Add the faces in the increasing order of neighbours
        for (label neiSearch = 0; neiSearch < nNeighbours; neiSearch++)
        {
            // Find the lowest neighbour which is still valid
            label nextNei = -1;
            label minNei = cellPolys_.size();
 
            forAll(neiCells, ncI)
            {
                if (neiCells[ncI] > -1 && neiCells[ncI] < minNei)
                {
                    nextNei = ncI;
                    minNei = neiCells[ncI];
                }
            }
 
            if (nextNei > -1)
            {
                // Add the face to the list of faces
                meshFaces_[nInternalFaces_] = curFaces[nextNei];
 
                // Mark for owner
                cellPolys_[celli][nextNei] = nInternalFaces_;
 
                // Mark for neighbour
                cellPolys_[neiCells[nextNei]][faceOfNeiCell[nextNei]] =
                    nInternalFaces_;
 
                // Stop the neighbour from being used again
                neiCells[nextNei] = -1;
 
                // Increment number of faces counter
                nInternalFaces_++;
            }
            else
            {
              FatalErrorInFunction
                  << "Error in internal face insertion"
                  << abort(FatalError);
            }
        }
    }
 
#ifdef DEBUG_CELLPOLY
    Info<< "cellPolys = " << cellPolys_ << endl;
#endif
 
    // don't reset the size of internal faces, because more faces will be
    // added in createPolyBoundary()
}
 
 
// ************************************************************************* //