foam-extend4.1-coherent-io/applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | foam-extend: Open Source CFD
   \\    /   O peration     | Version:     4.0
    \\  /    A nd           | Web:         http://www.foam-extend.org
     \\/     M anispulation  |
-------------------------------------------------------------------------------
License
    This file is part of foam-extend.

    foam-extend 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.

    foam-extend 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 foam-extend.  If not, see <http://www.gnu.org/licenses/>.

Application
    renumberMesh

Description
    Renumbers the cell list in order to reduce the bandwidth, reading and
    renumbering all fields from all the time directories.

\*---------------------------------------------------------------------------*/

#include "argList.H"
#include "IOobjectList.H"
#include "fvMesh.H"
#include "directTopoChange.H"
#include "ReadFields.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "bandCompression.H"
#include "faceSet.H"
#include "SortableList.H"
#include "decompositionMethod.H"
#include "fvMeshSubset.H"
#include "zeroGradientFvPatchFields.H"

using namespace Foam;


// Calculate band of matrix
label getBand(const labelList& owner, const labelList& neighbour)
{
    label band = 0;

    forAll(neighbour, faceI)
    {
        label diff = neighbour[faceI] - owner[faceI];

        if (diff > band)
        {
            band = diff;
        }
    }
    return band;
}


// Return new to old cell numbering
labelList regionBandCompression
(
    const fvMesh& mesh,
    const labelList& cellToRegion
)
{
    Pout<< "Determining cell order:" << endl;

    labelList cellOrder(cellToRegion.size());

    label nRegions = max(cellToRegion)+1;

    labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));

    label cellI = 0;

    forAll(regionToCells, regionI)
    {
        Pout<< "    region " << regionI << " starts at " << cellI << endl;

        // Per region do a reordering.
        fvMeshSubset meshSubset
        (
            IOobject
            (
                "set",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh
        );
        meshSubset.setLargeCellSubset(cellToRegion, regionI);
        const fvMesh& subMesh = meshSubset.subMesh();
        labelList subCellOrder(bandCompression(subMesh.cellCells()));

        const labelList& cellMap = meshSubset.cellMap();

        forAll(subCellOrder, i)
        {
            cellOrder[cellI++] = cellMap[subCellOrder[i]];
        }
    }
    Pout<< endl;

    return cellOrder;
}


// Determine face order such that inside region faces are sorted
// upper-triangular but inbetween region faces are handled like boundary faces.
labelList regionFaceOrder
(
    const primitiveMesh& mesh,
    const labelList& cellOrder,     // new to old cell
    const labelList& cellToRegion   // old cell to region
)
{
    Pout<< "Determining face order:" << endl;

    labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));

    labelList oldToNewFace(mesh.nFaces(), -1);

    label newFaceI = 0;

    label prevRegion = -1;

    forAll (cellOrder, newCellI)
    {
        label oldCellI = cellOrder[newCellI];

        if (cellToRegion[oldCellI] != prevRegion)
        {
            prevRegion = cellToRegion[oldCellI];
            Pout<< "    region " << prevRegion << " internal faces start at "
                << newFaceI << endl;
        }

        const cell& cFaces = mesh.cells()[oldCellI];

        SortableList<label> nbr(cFaces.size());

        forAll(cFaces, i)
        {
            label faceI = cFaces[i];

            if (mesh.isInternalFace(faceI))
            {
                // Internal face. Get cell on other side.
                label nbrCellI = reverseCellOrder[mesh.faceNeighbour()[faceI]];
                if (nbrCellI == newCellI)
                {
                    nbrCellI = reverseCellOrder[mesh.faceOwner()[faceI]];
                }

                if (cellToRegion[oldCellI] != cellToRegion[cellOrder[nbrCellI]])
                {
                    // Treat like external face. Do later.
                    nbr[i] = -1;
                }
                else if (newCellI < nbrCellI)
                {
                    // CellI is master
                    nbr[i] = nbrCellI;
                }
                else
                {
                    // nbrCell is master. Let it handle this face.
                    nbr[i] = -1;
                }
            }
            else
            {
                // External face. Do later.
                nbr[i] = -1;
            }
        }

        nbr.sort();

        forAll(nbr, i)
        {
            if (nbr[i] != -1)
            {
                oldToNewFace[cFaces[nbr.indices()[i]]] = newFaceI++;
            }
        }
    }

    // Do region interfaces
    label nRegions = max(cellToRegion)+1;
    {
        // Sort in increasing region
        SortableList<label> sortKey(mesh.nFaces(), labelMax);

        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
        {
            label ownRegion = cellToRegion[mesh.faceOwner()[faceI]];
            label neiRegion = cellToRegion[mesh.faceNeighbour()[faceI]];

            if (ownRegion != neiRegion)
            {
                sortKey[faceI] =
                    min(ownRegion, neiRegion)*nRegions
                   +max(ownRegion, neiRegion);
            }
        }
        sortKey.sort();

        // Extract.
        label prevKey = -1;
        forAll(sortKey, i)
        {
            label key = sortKey[i];

            if (key == labelMax)
            {
                break;
            }

            if (prevKey != key)
            {
                Pout<< "    faces inbetween region " << key/nRegions
                    << " and " << key%nRegions
                    << " start at " << newFaceI << endl;
                prevKey = key;
            }

            oldToNewFace[sortKey.indices()[i]] = newFaceI++;
        }
    }

    // Leave patch faces intact.
    for (label faceI = newFaceI; faceI < mesh.nFaces(); faceI++)
    {
        oldToNewFace[faceI] = faceI;
    }


    // Check done all faces.
    forAll(oldToNewFace, faceI)
    {
        if (oldToNewFace[faceI] == -1)
        {
            FatalErrorIn
            (
                "polyDualMesh::getFaceOrder"
                "(const labelList&, const labelList&, const label) const"
            )   << "Did not determine new position"
                << " for face " << faceI
                << abort(FatalError);
        }
    }
    Pout<< endl;

    return invert(mesh.nFaces(), oldToNewFace);
}


// cellOrder: old cell for every new cell
// faceOrder: old face for every new face. Ordering of boundary faces not
// changed.
autoPtr<mapPolyMesh> reorderMesh
(
    polyMesh& mesh,
    const labelList& cellOrder,
    const labelList& faceOrder
)
{
    labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
    labelList reverseFaceOrder(invert(faceOrder.size(), faceOrder));

    faceList newFaces(reorder(reverseFaceOrder, mesh.faces()));
    labelList newOwner
    (
        renumber
        (
            reverseCellOrder,
            reorder(reverseFaceOrder, mesh.faceOwner())
        )
    );
    labelList newNeighbour
    (
        renumber
        (
            reverseCellOrder,
            reorder(reverseFaceOrder, mesh.faceNeighbour())
        )
    );

    // Check if any faces need swapping.
    forAll(newNeighbour, faceI)
    {
        label own = newOwner[faceI];
        label nei = newNeighbour[faceI];

        if (nei < own)
        {
            newFaces[faceI] = newFaces[faceI].reverseFace();
            Swap(newOwner[faceI], newNeighbour[faceI]);
        }
    }

    const polyBoundaryMesh& patches = mesh.boundaryMesh();
    labelList patchSizes(patches.size());
    labelList patchStarts(patches.size());
    labelList oldPatchNMeshPoints(patches.size());
    labelListList patchPointMap(patches.size());
    forAll(patches, patchI)
    {
        patchSizes[patchI] = patches[patchI].size();
        patchStarts[patchI] = patches[patchI].start();
        oldPatchNMeshPoints[patchI] = patches[patchI].nPoints();
        patchPointMap[patchI] = identity(patches[patchI].nPoints());
    }

    mesh.resetPrimitives
    (
        Xfer<pointField>::null(),
        xferMove(newFaces),
        xferMove(newOwner),
        xferMove(newNeighbour),
        patchSizes,
        patchStarts
    );

    return autoPtr<mapPolyMesh>
    (
        new mapPolyMesh
        (
            mesh,                       //const polyMesh& mesh,
            mesh.nPoints(),             // nOldPoints,
            mesh.nFaces(),              // nOldFaces,
            mesh.nCells(),              // nOldCells,
            identity(mesh.nPoints()),   // pointMap,
            List<objectMap>(0),         // pointsFromPoints,
            faceOrder,                  // faceMap,
            List<objectMap>(0),         // facesFromPoints,
            List<objectMap>(0),         // facesFromEdges,
            List<objectMap>(0),         // facesFromFaces,
            cellOrder,                  // cellMap,
            List<objectMap>(0),         // cellsFromPoints,
            List<objectMap>(0),         // cellsFromEdges,
            List<objectMap>(0),         // cellsFromFaces,
            List<objectMap>(0),         // cellsFromCells,
            identity(mesh.nPoints()),   // reversePointMap,
            reverseFaceOrder,           // reverseFaceMap,
            reverseCellOrder,           // reverseCellMap,
            labelHashSet(0),            // flipFaceFlux,
            patchPointMap,              // patchPointMap,
            labelListList(0),           // pointZoneMap,
            labelListList(0),           // faceZonePointMap,
            labelListList(0),           // faceZoneFaceMap,
            labelListList(0),           // cellZoneMap,
            pointField(0),              // preMotionPoints,
            patchStarts,                // oldPatchStarts,
            oldPatchNMeshPoints         // oldPatchNMeshPoints
        )
    );
}


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

int main(int argc, char *argv[])
{
    argList::validOptions.insert("blockOrder", "");
    argList::validOptions.insert("writeMaps", "");
    argList::validOptions.insert("overwrite", "");

#   include "addTimeOptions.H"

#   include "setRootCase.H"
#   include "createTime.H"
    runTime.functionObjects().off();

    // Get times list
    instantList Times = runTime.times();

    // set startTime and endTime depending on -time and -latestTime options
#   include "checkTimeOptions.H"

    runTime.setTime(Times[startTime], startTime);

#   include "createMesh.H"
    const word oldInstance = mesh.pointsInstance();

    const bool blockOrder = args.optionFound("blockOrder");
    if (blockOrder)
    {
        Info<< "Ordering cells into regions (using decomposition);"
            << " ordering faces into region-internal and region-external." << nl
            << endl;
    }

    const bool writeMaps = args.optionFound("writeMaps");

    if (writeMaps)
    {
        Info<< "Writing renumber maps (new to old) to polyMesh." << nl
            << endl;
    }

    bool overwrite = args.optionFound("overwrite");

    label band = getBand(mesh.faceOwner(), mesh.faceNeighbour());

    Info<< "Mesh size: " << returnReduce(mesh.nCells(), sumOp<label>()) << nl
        << "Band before renumbering: "
        << returnReduce(band, maxOp<label>()) << nl << endl;

    // Read objects in time directory
    IOobjectList objects(mesh, runTime.timeName());

    // Read vol fields.

    PtrList<volScalarField> vsFlds;
    ReadFields(mesh, objects, vsFlds);

    PtrList<volVectorField> vvFlds;
    ReadFields(mesh, objects, vvFlds);

    PtrList<volSphericalTensorField> vstFlds;
    ReadFields(mesh, objects, vstFlds);

    PtrList<volSymmTensorField> vsymtFlds;
    ReadFields(mesh, objects, vsymtFlds);

    PtrList<volTensorField> vtFlds;
    ReadFields(mesh, objects, vtFlds);

    // Read surface fields.

    PtrList<surfaceScalarField> ssFlds;
    ReadFields(mesh, objects, ssFlds);

    PtrList<surfaceVectorField> svFlds;
    ReadFields(mesh, objects, svFlds);

    PtrList<surfaceSphericalTensorField> sstFlds;
    ReadFields(mesh, objects, sstFlds);

    PtrList<surfaceSymmTensorField> ssymtFlds;
    ReadFields(mesh, objects, ssymtFlds);

    PtrList<surfaceTensorField> stFlds;
    ReadFields(mesh, objects, stFlds);


    autoPtr<mapPolyMesh> map;

    if (blockOrder)
    {
        // Renumbering in two phases. Should be done in one so mapping of
        // fields is done correctly!

        // Read decomposePar dictionary
        IOdictionary decomposeDict
        (
            IOobject
            (
                "decomposeParDict",
                runTime.system(),
                mesh,
                IOobject::MUST_READ,
                IOobject::NO_WRITE
            )
        );
        autoPtr<decompositionMethod> decomposePtr = decompositionMethod::New
        (
            decomposeDict,
            mesh
        );

        labelList cellToRegion(decomposePtr().decompose(mesh.cellCentres()));

        // For debugging: write out region
        {
            volScalarField cellDist
            (
                IOobject
                (
                    "cellDist",
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE,
                    false
                ),
                mesh,
                dimensionedScalar("cellDist", dimless, 0),
                zeroGradientFvPatchScalarField::typeName
            );

            forAll(cellToRegion, cellI)
            {
               cellDist[cellI] = cellToRegion[cellI];
            }

            cellDist.write();

            Info<< nl << "Written decomposition as volScalarField to "
                << cellDist.name() << " for use in postprocessing."
                << nl << endl;
        }

        // Use block based renumbering.
        //labelList cellOrder(bandCompression(mesh.cellCells()));
        labelList cellOrder(regionBandCompression(mesh, cellToRegion));

        // Determine new to old face order with new cell numbering
        labelList faceOrder
        (
            regionFaceOrder
            (
                mesh,
                cellOrder,
                cellToRegion
            )
        );

        if (!overwrite)
        {
            runTime++;
        }

        // Change the mesh.
        map = reorderMesh(mesh, cellOrder, faceOrder);
    }
    else
    {
        // Use built-in renumbering.

        directTopoChange meshMod(mesh);

        if (!overwrite)
        {
            runTime++;
        }

        // Change the mesh.
        map = meshMod.changeMesh
        (
            mesh,
            false,      // inflate
            true,       // parallel sync
            true,       // cell ordering
            false       // point ordering
        );
    }

    // Update fields
    mesh.updateMesh(map);

    // Move mesh (since morphing might not do this)
    if (map().hasMotionPoints())
    {
        mesh.movePoints(map().preMotionPoints());
    }


    band = getBand(mesh.faceOwner(), mesh.faceNeighbour());

    Info<< "Band after renumbering: "
        << returnReduce(band, maxOp<label>()) << nl << endl;

    // Removed.  HJ, 23/Sep/2010
//     if (orderPoints)

    if (overwrite)
    {
        mesh.setInstance(oldInstance);
    }
    Info<< "Writing mesh to " << runTime.timeName() << endl;

    mesh.write();

    if (writeMaps)
    {
        labelIOList
        (
            IOobject
            (
                "cellMap",
                mesh.facesInstance(),
                polyMesh::meshSubDir,
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            map().cellMap()
        ).write();

        labelIOList
        (
            IOobject
            (
                "faceMap",
                mesh.facesInstance(),
                polyMesh::meshSubDir,
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            map().faceMap()
        ).write();

        labelIOList
        (
            IOobject
            (
                "pointMap",
                mesh.facesInstance(),
                polyMesh::meshSubDir,
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            map().pointMap()
        ).write();
    }

    Info<< "\nEnd.\n" << endl;

    return 0;
}


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