viewFactorsGen.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-2020 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/>.
 
Application
    viewFactorsGen
 
Group
    grpPreProcessingUtilities
 
Description
    View factors are calculated based on a face agglomeration array
    (finalAgglom generated by faceAgglomerate utility).
 
    Each view factor between the agglomerated faces i and j (Fij) is calculated
    using a double integral of the sub-areas composing the agglomeration.
 
    The patches involved in the view factor calculation are taken from the
    boundary file and should be part on the group viewFactorWall. ie.:
 
    floor
    {
        type            wall;
        inGroups        2(wall viewFactorWall);
        nFaces          100;
        startFace       3100;
    }
 
\*---------------------------------------------------------------------------*/
 
#include "argList.H"
#include "fvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "distributedTriSurfaceMesh.H"
#include "meshTools.H"
 
#include "uindirectPrimitivePatch.H"
#include "DynamicField.H"
#include "unitConversion.H"
 
#include "scalarMatrices.H"
#include "labelListIOList.H"
#include "scalarListIOList.H"
 
#include "singleCellFvMesh.H"
 
#include "IOmapDistribute.H"
 
using namespace Foam;
 
 
triSurface triangulate
(
    const polyBoundaryMesh& bMesh,
    const labelHashSet& includePatches,
    const labelListIOList& finalAgglom,
    labelList& triSurfaceToAgglom,
    const globalIndex& globalNumbering,
    const polyBoundaryMesh& coarsePatches
)
{
    const polyMesh& mesh = bMesh.mesh();
 
    // Storage for surfaceMesh. Size estimate.
    DynamicList<labelledTri> triangles(mesh.nBoundaryFaces());
 
    label newPatchI = 0;
    label localTriFaceI = 0;
 
    for (const label patchI : includePatches)
    {
        const polyPatch& patch = bMesh[patchI];
        const pointField& points = patch.points();
 
        label nTriTotal = 0;
 
        forAll(patch, patchFaceI)
        {
            const face& f = patch[patchFaceI];
 
            faceList triFaces(f.nTriangles(points));
 
            label nTri = 0;
 
            f.triangles(points, nTri, triFaces);
 
            forAll(triFaces, triFaceI)
            {
                const face& f = triFaces[triFaceI];
 
                triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
 
                nTriTotal++;
 
                triSurfaceToAgglom[localTriFaceI++] = globalNumbering.toGlobal
                (
                    Pstream::myProcNo(),
                    finalAgglom[patchI][patchFaceI]
                  + coarsePatches[patchI].start()
                );
            }
        }
 
        newPatchI++;
    }
 
    //striSurfaceToAgglom.resize(localTriFaceI-1);
 
    triangles.shrink();
 
    // Create globally numbered tri surface
    triSurface rawSurface(triangles, mesh.points());
 
    // Create locally numbered tri surface
    triSurface surface
    (
        rawSurface.localFaces(),
        rawSurface.localPoints()
    );
 
    // Add patch names to surface
    surface.patches().setSize(newPatchI);
 
    newPatchI = 0;
 
    for (const label patchI : includePatches)
    {
        const polyPatch& patch = bMesh[patchI];
 
        surface.patches()[newPatchI].index() = patchI;
        surface.patches()[newPatchI].name() = patch.name();
        surface.patches()[newPatchI].geometricType() = patch.type();
 
        newPatchI++;
    }
 
    return surface;
}
 
 
void writeRays
(
    const fileName& fName,
    const pointField& compactCf,
    const pointField& myFc,
    const labelListList& visibleFaceFaces
)
{
    OFstream str(fName);
    label vertI = 0;
 
    Pout<< "Dumping rays to " << str.name() << endl;
 
    forAll(myFc, faceI)
    {
        const labelList visFaces = visibleFaceFaces[faceI];
        forAll(visFaces, faceRemote)
        {
            label compactI = visFaces[faceRemote];
            const point& remoteFc = compactCf[compactI];
 
            meshTools::writeOBJ(str, myFc[faceI]);
            vertI++;
            meshTools::writeOBJ(str, remoteFc);
            vertI++;
            str << "l " << vertI-1 << ' ' << vertI << nl;
        }
    }
    str.flush();
}
 
 
scalar calculateViewFactorFij
(
    const vector& i,
    const vector& j,
    const vector& dAi,
    const vector& dAj
)
{
    vector r = i - j;
    scalar rMag = mag(r);
 
    if (rMag > SMALL)
    {
        scalar dAiMag = mag(dAi);
        scalar dAjMag = mag(dAj);
 
        vector ni = dAi/dAiMag;
        vector nj = dAj/dAjMag;
        scalar cosThetaJ = mag(nj & r)/rMag;
        scalar cosThetaI = mag(ni & r)/rMag;
 
        return
        (
            (cosThetaI*cosThetaJ*dAjMag*dAiMag)
           /(sqr(rMag)*constant::mathematical::pi)
        );
    }
    else
    {
        return 0;
    }
}
 
 
void insertMatrixElements
(
    const globalIndex& globalNumbering,
    const label fromProcI,
    const labelListList& globalFaceFaces,
    const scalarListList& viewFactors,
    scalarSquareMatrix& matrix
)
{
    forAll(viewFactors, faceI)
    {
        const scalarList& vf = viewFactors[faceI];
        const labelList& globalFaces = globalFaceFaces[faceI];
 
        label globalI = globalNumbering.toGlobal(fromProcI, faceI);
        forAll(globalFaces, i)
        {
            matrix[globalI][globalFaces[i]] = vf[i];
        }
    }
}
 
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Calculate view factors from face agglomeration array."
        " The finalAgglom generated by faceAgglomerate utility."
    );
 
    #include "addRegionOption.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createNamedMesh.H"
 
    // Read view factor dictionary
    IOdictionary viewFactorDict
    (
       IOobject
       (
            "viewFactorsDict",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
       )
    );
 
    const word viewFactorWall("viewFactorWall");
 
    const bool writeViewFactors =
        viewFactorDict.getOrDefault("writeViewFactorMatrix", false);
 
    const bool dumpRays =
        viewFactorDict.getOrDefault("dumpRays", false);
 
    const label debug = viewFactorDict.getOrDefault<label>("debug", 0);
 
    // Read agglomeration map
    labelListIOList finalAgglom
    (
        IOobject
        (
            "finalAgglom",
            mesh.facesInstance(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE,
            false
        )
    );
 
    // Create the coarse mesh  using agglomeration
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    if (debug)
    {
        Pout << "\nCreating single cell mesh..." << endl;
    }
 
    singleCellFvMesh coarseMesh
    (
        IOobject
        (
            "coarse:" + mesh.name(),
            runTime.timeName(),
            runTime,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        mesh,
        finalAgglom
    );
 
    if (debug)
    {
        Pout << "\nCreated single cell mesh..." << endl;
    }
 
 
    // Calculate total number of fine and coarse faces
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    label nCoarseFaces = 0;      //total number of coarse faces
    label nFineFaces = 0;        //total number of fine faces
 
    const polyBoundaryMesh& patches = mesh.boundaryMesh();
    const polyBoundaryMesh& coarsePatches = coarseMesh.boundaryMesh();
 
    labelList viewFactorsPatches(patches.indices(viewFactorWall));
    for (const label patchi : viewFactorsPatches)
    {
        nCoarseFaces += coarsePatches[patchi].size();
        nFineFaces += patches[patchi].size();
    }
 
    // total number of coarse faces
    label totalNCoarseFaces = nCoarseFaces;
 
    reduce(totalNCoarseFaces, sumOp<label>());
 
    if (Pstream::master())
    {
        Info << "\nTotal number of coarse faces: "<< totalNCoarseFaces << endl;
    }
 
    if (Pstream::master() && debug)
    {
        Pout << "\nView factor patches included in the calculation : "
             << viewFactorsPatches << endl;
    }
 
    // Collect local Cf and Sf on coarse mesh
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    DynamicList<point> localCoarseCf(nCoarseFaces);
    DynamicList<point> localCoarseSf(nCoarseFaces);
    DynamicList<label> localAgg(nCoarseFaces);
    labelHashSet includePatches;
 
    for (const label patchID : viewFactorsPatches)
    {
        const polyPatch& pp = patches[patchID];
        const labelList& agglom = finalAgglom[patchID];
 
        includePatches.insert(patchID);
 
        if (agglom.size() > 0)
        {
            label nAgglom = max(agglom)+1;
            labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
            const labelList& coarsePatchFace =
                coarseMesh.patchFaceMap()[patchID];
 
            const pointField& coarseCf =
                coarseMesh.Cf().boundaryField()[patchID];
            const pointField& coarseSf =
                coarseMesh.Sf().boundaryField()[patchID];
 
            forAll(coarseCf, faceI)
            {
                point cf = coarseCf[faceI];
 
                const label coarseFaceI = coarsePatchFace[faceI];
                const labelList& fineFaces = coarseToFine[coarseFaceI];
                const label agglomI =
                    agglom[fineFaces[0]] + coarsePatches[patchID].start();
 
                // Construct single face
                uindirectPrimitivePatch upp
                (
                    UIndirectList<face>(pp, fineFaces),
                    pp.points()
                );
 
                List<point> availablePoints
                (
                    upp.faceCentres().size()
                + upp.localPoints().size()
                );
 
                SubList<point>
                (
                    availablePoints,
                    upp.faceCentres().size()
                ) = upp.faceCentres();
 
                SubList<point>
                (
                    availablePoints,
                    upp.localPoints().size(),
                    upp.faceCentres().size()
                ) = upp.localPoints();
 
                point cfo = cf;
                scalar dist = GREAT;
                forAll(availablePoints, iPoint)
                {
                    point cfFine = availablePoints[iPoint];
                    if (mag(cfFine-cfo) < dist)
                    {
                        dist = mag(cfFine-cfo);
                        cf = cfFine;
                    }
                }
 
                point sf = coarseSf[faceI];
                localCoarseCf.append(cf);
                localCoarseSf.append(sf);
                localAgg.append(agglomI);
 
            }
        }
    }
 
    // Distribute local coarse Cf and Sf for shooting rays
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    List<pointField> remoteCoarseCf(Pstream::nProcs());
    List<pointField> remoteCoarseSf(Pstream::nProcs());
    List<labelField> remoteCoarseAgg(Pstream::nProcs());
 
    remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
    remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
    remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
 
    Pstream::gatherList(remoteCoarseCf);
    Pstream::scatterList(remoteCoarseCf);
    Pstream::gatherList(remoteCoarseSf);
    Pstream::scatterList(remoteCoarseSf);
    Pstream::gatherList(remoteCoarseAgg);
    Pstream::scatterList(remoteCoarseAgg);
 
 
    globalIndex globalNumbering(nCoarseFaces);
 
    // Set up searching engine for obstacles
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    #include "searchingEngine.H"
 
    // Determine rays between coarse face centres
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    DynamicList<label> rayStartFace(nCoarseFaces + 0.01*nCoarseFaces);
 
    DynamicList<label> rayEndFace(rayStartFace.size());
 
 
    // Return rayStartFace in local index and rayEndFace in global index
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    #include "shootRays.H"
 
    // Calculate number of visible faces from local index
    labelList nVisibleFaceFaces(nCoarseFaces, Zero);
 
    forAll(rayStartFace, i)
    {
        nVisibleFaceFaces[rayStartFace[i]]++;
    }
 
    labelListList visibleFaceFaces(nCoarseFaces);
 
    label nViewFactors = 0;
    forAll(nVisibleFaceFaces, faceI)
    {
        visibleFaceFaces[faceI].setSize(nVisibleFaceFaces[faceI]);
        nViewFactors += nVisibleFaceFaces[faceI];
    }
 
    // - Construct compact numbering
    // - return map from remote to compact indices
    //   (per processor (!= myProcNo) a map from remote index to compact index)
    // - construct distribute map
    // - renumber rayEndFace into compact addressing
 
    List<Map<label>> compactMap(Pstream::nProcs());
 
    mapDistribute map(globalNumbering, rayEndFace, compactMap);
 
    // visibleFaceFaces has:
    //    (local face, local viewed face) = compact viewed face
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    nVisibleFaceFaces = 0;
    forAll(rayStartFace, i)
    {
        label faceI = rayStartFace[i];
        label compactI = rayEndFace[i];
        visibleFaceFaces[faceI][nVisibleFaceFaces[faceI]++] = compactI;
    }
 
 
    // Construct data in compact addressing
    // I need coarse Sf (Ai), fine Sf (dAi) and fine Cf(r) to calculate Fij
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    pointField compactCoarseCf(map.constructSize(), Zero);
    pointField compactCoarseSf(map.constructSize(), Zero);
    List<List<point>> compactFineSf(map.constructSize());
    List<List<point>> compactFineCf(map.constructSize());
 
    DynamicList<label> compactPatchId(map.constructSize());
 
    // Insert my coarse local values
    SubList<point>(compactCoarseSf, nCoarseFaces) = localCoarseSf;
    SubList<point>(compactCoarseCf, nCoarseFaces) = localCoarseCf;
 
    // Insert my fine local values
    label compactI = 0;
    forAll(viewFactorsPatches, i)
    {
        label patchID = viewFactorsPatches[i];
 
        const labelList& agglom = finalAgglom[patchID];
        if (agglom.size() > 0)
        {
            label nAgglom = max(agglom)+1;
            labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
            const labelList& coarsePatchFace =
                coarseMesh.patchFaceMap()[patchID];
 
            forAll(coarseToFine, coarseI)
            {
                compactPatchId.append(patchID);
                List<point>& fineCf = compactFineCf[compactI];
                List<point>& fineSf = compactFineSf[compactI++];
 
                const label coarseFaceI = coarsePatchFace[coarseI];
                const labelList& fineFaces = coarseToFine[coarseFaceI];
 
                fineCf.setSize(fineFaces.size());
                fineSf.setSize(fineFaces.size());
 
                fineCf = UIndirectList<point>
                (
                    mesh.Cf().boundaryField()[patchID],
                    coarseToFine[coarseFaceI]
                );
                fineSf = UIndirectList<point>
                (
                    mesh.Sf().boundaryField()[patchID],
                    coarseToFine[coarseFaceI]
                );
            }
        }
    }
 
    // Do all swapping
    map.distribute(compactCoarseSf);
    map.distribute(compactCoarseCf);
    map.distribute(compactFineCf);
    map.distribute(compactFineSf);
 
    map.distribute(compactPatchId);
 
 
    // Plot all rays between visible faces.
    if (dumpRays)
    {
        writeRays
        (
            runTime.path()/"allVisibleFaces.obj",
            compactCoarseCf,
            remoteCoarseCf[Pstream::myProcNo()],
            visibleFaceFaces
        );
    }
 
 
    // Fill local view factor matrix
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
    scalarListIOList F
    (
        IOobject
        (
            "F",
            mesh.facesInstance(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        nCoarseFaces
    );
 
    label totalPatches = coarsePatches.size();
    reduce(totalPatches, maxOp<label>());
 
    // Matrix sum in j(Fij) for each i (if enclosure sum = 1)
    scalarSquareMatrix sumViewFactorPatch
    (
        totalPatches,
        0.0
    );
 
    scalarList patchArea(totalPatches, Zero);
 
    if (Pstream::master())
    {
        Info<< "\nCalculating view factors..." << endl;
    }
 
    if (mesh.nSolutionD() == 3)
    {
        forAll(localCoarseSf, coarseFaceI)
        {
            const List<point>& localFineSf = compactFineSf[coarseFaceI];
            const vector Ai = sum(localFineSf);
            const List<point>& localFineCf = compactFineCf[coarseFaceI];
            const label fromPatchId = compactPatchId[coarseFaceI];
            patchArea[fromPatchId] += mag(Ai);
 
            const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
 
            forAll(visCoarseFaces, visCoarseFaceI)
            {
                F[coarseFaceI].setSize(visCoarseFaces.size());
                label compactJ = visCoarseFaces[visCoarseFaceI];
                const List<point>& remoteFineSj = compactFineSf[compactJ];
                const List<point>& remoteFineCj = compactFineCf[compactJ];
 
                const label toPatchId = compactPatchId[compactJ];
 
                scalar Fij = 0;
                forAll(localFineSf, i)
                {
                    const vector& dAi = localFineSf[i];
                    const vector& dCi = localFineCf[i];
 
                    forAll(remoteFineSj, j)
                    {
                        const vector& dAj = remoteFineSj[j];
                        const vector& dCj = remoteFineCj[j];
 
                        scalar dIntFij = calculateViewFactorFij
                        (
                            dCi,
                            dCj,
                            dAi,
                            dAj
                        );
 
                        Fij += dIntFij;
                    }
                }
                F[coarseFaceI][visCoarseFaceI] = Fij/mag(Ai);
                sumViewFactorPatch[fromPatchId][toPatchId] += Fij;
            }
        }
    }
    else if (mesh.nSolutionD() == 2)
    {
        const boundBox& box = mesh.bounds();
        const Vector<label>& dirs = mesh.geometricD();
        vector emptyDir = Zero;
        forAll(dirs, i)
        {
            if (dirs[i] == -1)
            {
                emptyDir[i] = 1.0;
            }
        }
 
        scalar wideBy2 = (box.span() & emptyDir)*2.0;
 
        forAll(localCoarseSf, coarseFaceI)
        {
            const vector& Ai = localCoarseSf[coarseFaceI];
            const vector& Ci = localCoarseCf[coarseFaceI];
            vector Ain = Ai/mag(Ai);
            vector R1i = Ci + (mag(Ai)/wideBy2)*(Ain ^ emptyDir);
            vector R2i = Ci - (mag(Ai)/wideBy2)*(Ain ^ emptyDir) ;
 
            const label fromPatchId = compactPatchId[coarseFaceI];
            patchArea[fromPatchId] += mag(Ai);
 
            const labelList& visCoarseFaces = visibleFaceFaces[coarseFaceI];
            forAll(visCoarseFaces, visCoarseFaceI)
            {
                F[coarseFaceI].setSize(visCoarseFaces.size());
                label compactJ = visCoarseFaces[visCoarseFaceI];
                const vector& Aj = compactCoarseSf[compactJ];
                const vector& Cj = compactCoarseCf[compactJ];
 
                const label toPatchId = compactPatchId[compactJ];
 
                vector Ajn = Aj/mag(Aj);
                vector R1j = Cj + (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
                vector R2j = Cj - (mag(Aj)/wideBy2)*(Ajn ^ emptyDir);
 
                scalar d1 = mag(R1i - R2j);
                scalar d2 = mag(R2i - R1j);
                scalar s1 = mag(R1i - R1j);
                scalar s2 = mag(R2i - R2j);
 
                scalar Fij = mag((d1 + d2) - (s1 + s2))/(4.0*mag(Ai)/wideBy2);
 
                F[coarseFaceI][visCoarseFaceI] = Fij;
                sumViewFactorPatch[fromPatchId][toPatchId] += Fij*mag(Ai);
            }
        }
    }
 
    if (Pstream::master())
    {
        Info << "Writing view factor matrix..." << endl;
    }
 
    // Write view factors matrix in listlist form
    F.write();
 
    reduce(sumViewFactorPatch, sumOp<scalarSquareMatrix>());
    reduce(patchArea, sumOp<scalarList>());
 
 
    if (Pstream::master() && debug)
    {
        forAll(viewFactorsPatches, i)
        {
            label patchI =  viewFactorsPatches[i];
            forAll(viewFactorsPatches, i)
            {
                label patchJ =  viewFactorsPatches[i];
                Info << "F" << patchI << patchJ << ": "
                     << sumViewFactorPatch[patchI][patchJ]/patchArea[patchI]
                     << endl;
            }
        }
    }
 
 
    if (writeViewFactors)
    {
        volScalarField viewFactorField
        (
            IOobject
            (
                "viewFactorField",
                mesh.time().timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar(dimless, Zero)
        );
 
        label compactI = 0;
 
        volScalarField::Boundary& vfbf = viewFactorField.boundaryFieldRef();
        forAll(viewFactorsPatches, i)
        {
            label patchID = viewFactorsPatches[i];
            const labelList& agglom = finalAgglom[patchID];
            if (agglom.size() > 0)
            {
                label nAgglom = max(agglom)+1;
                labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
                const labelList& coarsePatchFace =
                    coarseMesh.patchFaceMap()[patchID];
 
                forAll(coarseToFine, coarseI)
                {
                    const scalar Fij = sum(F[compactI]);
                    const label coarseFaceID = coarsePatchFace[coarseI];
                    const labelList& fineFaces = coarseToFine[coarseFaceID];
                    forAll(fineFaces, fineId)
                    {
                        const label faceID = fineFaces[fineId];
                        vfbf[patchID][faceID] = Fij;
                    }
                    compactI++;
                }
            }
        }
        viewFactorField.write();
    }
 
 
    // Invert compactMap (from processor+localface to compact) to go
    // from compact to processor+localface (expressed as a globalIndex)
    // globalIndex globalCoarFaceNum(coarseMesh.nFaces());
    labelList compactToGlobal(map.constructSize());
 
    // Local indices first (note: are not in compactMap)
    for (label i = 0; i < globalNumbering.localSize(); i++)
    {
        compactToGlobal[i] = globalNumbering.toGlobal(i);
    }
 
 
    forAll(compactMap, procI)
    {
        const Map<label>& localToCompactMap = compactMap[procI];
 
        forAllConstIters(localToCompactMap, iter)
        {
            compactToGlobal[*iter] = globalNumbering.toGlobal
            (
                procI,
                iter.key()
            );
        }
    }
 
 
    labelListList globalFaceFaces(visibleFaceFaces.size());
 
    // Create globalFaceFaces needed to insert view factors
    // in F to the global matrix Fmatrix
    forAll(globalFaceFaces, faceI)
    {
        globalFaceFaces[faceI] = renumber
        (
            compactToGlobal,
            visibleFaceFaces[faceI]
        );
    }
 
    labelListIOList IOglobalFaceFaces
    (
        IOobject
        (
            "globalFaceFaces",
            mesh.facesInstance(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        std::move(globalFaceFaces)
    );
    IOglobalFaceFaces.write();
 
 
    labelListIOList IOvisibleFaceFaces
    (
        IOobject
        (
            "visibleFaceFaces",
            mesh.facesInstance(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        std::move(visibleFaceFaces)
    );
    IOvisibleFaceFaces.write();
 
    IOmapDistribute IOmapDist
    (
        IOobject
        (
            "mapDist",
            mesh.facesInstance(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        std::move(map)
    );
 
    IOmapDist.write();
 
    Info<< "End\n" << endl;
    return 0;
}
 
 
// ************************************************************************* //