bandCompression.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-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 "bandCompression.H"
#include "bitSet.H"
#include "polyMesh.H"
#include "globalMeshData.H"
#include "CircularBuffer.H"
#include "CompactListList.H"
#include "DynamicList.H"
#include "IOstreams.H"
 
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
namespace
{
 
// Process connections with the Cuthill-McKee algorithm.
// The connections are CompactListList<label> or a labelListList.
template<class ConnectionListListType>
Foam::labelList cuthill_mckee_algorithm
(
    const ConnectionListListType& cellCellAddressing
)
{
    using namespace Foam;
 
    const label nOldCells(cellCellAddressing.size());
 
    // Which cells are visited/unvisited
    bitSet unvisited(nOldCells, true);
 
    // The new output order
    labelList newOrder(nOldCells);
 
 
    // Various work arrays
    // ~~~~~~~~~~~~~~~~~~~
 
    // Neighbour cells
    DynamicList<label> nbrCells;
 
    // Neighbour ordering
    DynamicList<label> nbrOrder;
 
    // Corresponding weights for neighbour cells
    DynamicList<label> weights;
 
    // FIFO buffer for renumbering.
    CircularBuffer<label> queuedCells(1024);
 
    label cellInOrder = 0;
 
    while (true)
    {
        // For a disconnected region find the lowest connected cell.
        label currCelli = -1;
        label minCount = labelMax;
 
        for (const label celli : unvisited)
        {
            const label nbrCount = cellCellAddressing[celli].size();
 
            if (minCount > nbrCount)
            {
                minCount = nbrCount;
                currCelli = celli;
            }
        }
 
        if (currCelli == -1)
        {
            break;
        }
 
 
        // Starting from currCelli - walk breadth-first
 
        queuedCells.push_back(currCelli);
 
        // Loop through queuedCells list. Add the first cell into the
        // cell order if it has not already been visited and ask for its
        // neighbours. If the neighbour in question has not been visited,
        // add it to the end of the queuedCells list
 
        while (!queuedCells.empty())
        {
            // Process as FIFO
            currCelli = queuedCells.front();
            queuedCells.pop_front();
 
            if (unvisited.test(currCelli))
            {
                // First visit...
                unvisited.unset(currCelli);
 
                // Add into cellOrder
                newOrder[cellInOrder] = currCelli;
                ++cellInOrder;
 
                // Add in increasing order of connectivity
 
                // 1. Count neighbours of unvisited neighbours
                nbrCells.clear();
                weights.clear();
 
                // Find if the neighbours have been visited
                const auto& neighbours = cellCellAddressing[currCelli];
 
                for (const label nbr : neighbours)
                {
                    const label nbrCount = cellCellAddressing[nbr].size();
 
                    if (unvisited.test(nbr))
                    {
                        // Not visited (or removed), add to the list
                        nbrCells.push_back(nbr);
                        weights.push_back(nbrCount);
                    }
                }
 
                // Resize DynamicList prior to sortedOrder
                nbrOrder.resize_nocopy(weights.size());
 
                // 2. Ascending order
                Foam::sortedOrder(weights, nbrOrder);
 
                // 3. Add to FIFO in sorted order
                for (const label nbrIdx : nbrOrder)
                {
                    queuedCells.push_back(nbrCells[nbrIdx]);
                }
            }
        }
    }
 
    // Debug:
    // - the peak capacity of queuedCells approximates the
    //   maximum intermediate bandwidth
    #if 0
    Pout<< "bandCompression: peak-capacity=" << queuedCells.capacity() << nl;
    #endif
 
    // Now we have new-to-old in newOrder.
    return newOrder;
}
 
} // End anonymous namespace
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
Foam::labelList Foam::meshTools::bandCompression
(
    const labelUList& cellCells,
    const labelUList& offsets
)
{
    // Protect against zero-sized offset list
    const label nOldCells = Foam::max(0, (offsets.size()-1));
 
    // Count number of neighbours
    labelList numNbrs(nOldCells, Foam::zero{});
    for (label celli = 0; celli < nOldCells; ++celli)
    {
        const label beg = offsets[celli];
        const label end = offsets[celli+1];
 
        for (label idx = beg; idx < end; ++idx)
        {
            ++numNbrs[celli];
            ++numNbrs[cellCells[idx]];
        }
    }
 
 
    // Which cells are visited/unvisited
    bitSet unvisited(nOldCells, true);
 
    // The new output order
    labelList newOrder(nOldCells);
 
 
    // Various work arrays
    // ~~~~~~~~~~~~~~~~~~~
 
    // Neighbour cells
    DynamicList<label> nbrCells;
 
    // Neighbour ordering
    DynamicList<label> nbrOrder;
 
    // Corresponding weights for neighbour cells
    DynamicList<label> weights;
 
    // FIFO buffer for renumbering.
    CircularBuffer<label> queuedCells(1024);
 
 
    label cellInOrder = 0;
 
    while (true)
    {
        // Find lowest connected cell that has not been visited yet
        label currCelli = -1;
        label minCount = labelMax;
 
        for (const label celli : unvisited)
        {
            const label nbrCount = numNbrs[celli];
 
            if (minCount > nbrCount)
            {
                minCount = nbrCount;
                currCelli = celli;
            }
        }
 
        if (currCelli == -1)
        {
            break;
        }
 
 
        // Starting from currCellii - walk breadth-first
 
        queuedCells.push_back(currCelli);
 
        // loop through the nextCell list. Add the first cell into the
        // cell order if it has not already been visited and ask for its
        // neighbours. If the neighbour in question has not been visited,
        // add it to the end of the nextCell list
 
        // Loop through queuedCells list. Add the first cell into the
        // cell order if it has not already been visited and ask for its
        // neighbours. If the neighbour in question has not been visited,
        // add it to the end of the queuedCells list
 
        while (!queuedCells.empty())
        {
            // Process as FIFO
            currCelli = queuedCells.front();
            queuedCells.pop_front();
 
            if (unvisited.test(currCelli))
            {
                // First visit...
                unvisited.unset(currCelli);
 
                // Add into cellOrder
                newOrder[cellInOrder] = currCelli;
                ++cellInOrder;
 
                // Add in increasing order of connectivity
 
                // 1. Count neighbours of unvisited neighbours
                nbrCells.clear();
                weights.clear();
 
                const label beg = offsets[currCelli];
                const label end = offsets[currCelli+1];
 
                for (label idx = beg; idx < end; ++idx)
                {
                    const label nbr = cellCells[idx];
                    const label nbrCount = numNbrs[nbr];
 
                    if (unvisited.test(nbr))
                    {
                        // Not visited (or removed), add to the list
                        nbrCells.push_back(nbr);
                        weights.push_back(nbrCount);
                    }
                }
 
                // Resize DynamicList prior to sortedOrder
                nbrOrder.resize_nocopy(weights.size());
 
                // 2. Ascending order
                Foam::sortedOrder(weights, nbrOrder);
 
                // 3. Add to FIFO in sorted order
                for (const label nbrIdx : nbrOrder)
                {
                    queuedCells.push_back(nbrCells[nbrIdx]);
                }
            }
        }
    }
 
    // Debug:
    // - the peak capacity of queuedCells approximates the
    //   maximum intermediate bandwidth
    #if 0
    Pout<< "bandCompression: peak-capacity=" << queuedCells.capacity() << nl;
    #endif
 
    // Now we have new-to-old in newOrder.
    return newOrder;
}
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
Foam::labelList Foam::meshTools::bandCompression(const polyMesh& mesh)
{
    // Local mesh connectivity
    CompactListList<label> cellCells;
    globalMeshData::calcCellCells(mesh, cellCells);
 
    return cuthill_mckee_algorithm(cellCells);
}
 
 
Foam::labelList Foam::meshTools::bandCompression
(
    const CompactListList<label>& cellCellAddressing
)
{
    return cuthill_mckee_algorithm(cellCellAddressing);
}
 
 
Foam::labelList Foam::meshTools::bandCompression
(
    const labelListList& cellCellAddressing
)
{
    return cuthill_mckee_algorithm(cellCellAddressing);
}
 
 
// ************************************************************************* //