/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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

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 agglomaration.

    The patches involved in the view factor calculation are taken from the qr
    volScalarField (radiative flux) when is greyDiffusiveRadiationViewFactor
    otherwise they are not included.

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

#include "argList.H"
#include "Time.H"
#include "singleCellFvMesh.H"
#include "volFields.H"
#include "surfaceFields.H"
#include "fixedValueFvPatchFields.H"
#include "distributedTriSurfaceMesh.H"
#include "symmetryPolyPatch.H"
#include "symmetryPlanePolyPatch.H"
#include "wedgePolyPatch.H"
#include "meshTools.H"
#include "uindirectPrimitivePatch.H"
#include "DynamicField.H"
#include "scalarListIOList.H"
#include "polygonTriangulate.H"
#include "vtkWritePolyData.H"
#include "scalarMatrices.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.nFaces() - mesh.nInternalFaces()
    );

    label newPatchi = 0;
    label localTriFacei = 0;

    polygonTriangulate triEngine;

    forAllConstIter(labelHashSet, includePatches, iter)
    {
        const label patchi = iter.key();
        const polyPatch& patch = bMesh[patchi];
        const pointField& points = patch.points();

        forAll(patch, patchFacei)
        {
            const face& f = patch[patchFacei];

            triEngine.triangulate(UIndirectList<point>(points, f));

            forAll(triEngine.triPoints(), triFacei)
            {
                triangles.append
                (
                    labelledTri
                    (
                        triEngine.triPoints(triFacei, f),
                        newPatchi
                    )
                );

                triSurfaceToAgglom[localTriFacei++] = globalNumbering.toGlobal
                (
                    Pstream::myProcNo(),
                    finalAgglom[patchi][patchFacei]
                  + coarsePatches[patchi].start()
                );
            }
        }

        newPatchi++;
    }

    triSurfaceToAgglom.resize(localTriFacei);

    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;

    forAllConstIter(labelHashSet, includePatches, iter)
    {
        const label patchi = iter.key();
        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
)
{
    DynamicList<point> allPoints;
    allPoints.append(myFc);
    allPoints.append(compactCf);

    DynamicList<labelPair> rays;
    forAll(myFc, facei)
    {
        const labelList visFaces = visibleFaceFaces[facei];
        forAll(visFaces, faceRemote)
        {
            rays.append
            (
                labelPair(facei, myFc.size() + visFaces[faceRemote])
            );
        }
    }

    Pout<< "\nDumping rays to " << fName + ".vtk" << endl;

    vtkWritePolyData::write
    (
        fName + ".vtk",
        fName.name(),
        false,
        allPoints,
        labelList(),
        rays,
        faceList()
    );
}


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[])
{
    #include "addRegionOption.H"
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createNamedMesh.H"

    const polyBoundaryMesh& patches = mesh.boundaryMesh();

    forAll(patches, patchi)
    {
        const polyPatch& pp = patches[patchi];
        if
        (
            isA<cyclicTransform>(pp)
         || isA<symmetryPolyPatch>(pp)
         || isA<symmetryPlanePolyPatch>(pp)
         || isA<wedgePolyPatch>(pp)
        )
        {
            FatalErrorIn(args.executable()) << args.executable()
                << " does not currently support transforming patches: "
                   "cyclic, symmetry and wedge."
                << exit(FatalError);
        }
    }

    // Read view factor dictionary
    IOdictionary viewFactorDict
    (
       IOobject
       (
            "viewFactorsDict",
            runTime.system(),
            mesh,
            IOobject::MUST_READ_IF_MODIFIED,
            IOobject::NO_WRITE
       )
    );

    const bool writeViewFactors =
        viewFactorDict.lookupOrDefault<bool>("writeViewFactorMatrix", false);

    const bool dumpRays =
        viewFactorDict.lookupOrDefault<bool>("dumpRays", false);

    const label debug = viewFactorDict.lookupOrDefault<label>("debug", 0);

    volScalarField qr
    (
        IOobject
        (
            "qr",
            runTime.name(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        ),
        mesh
    );

    // 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.name(),
            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& coarsePatches = coarseMesh.boundaryMesh();

    labelList viewFactorsPatches(patches.size());

    const volScalarField::Boundary& qrb = qr.boundaryField();

    label count = 0;
    forAll(qrb, patchi)
    {
        const polyPatch& pp = patches[patchi];
        const fvPatchScalarField& qrpI = qrb[patchi];

        if ((isA<fixedValueFvPatchScalarField>(qrpI)) && (pp.size() > 0))
        {
            viewFactorsPatches[count] = qrpI.patch().index();
            nCoarseFaces += coarsePatches[patchi].size();
            nFineFaces += patches[patchi].size();
            count ++;
        }
    }

    viewFactorsPatches.resize(count);

    // 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);

    forAll(viewFactorsPatches, i)
    {
        const label patchID = viewFactorsPatches[i];

        const polyPatch& pp = patches[patchID];
        const labelList& agglom = finalAgglom[patchID];
        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];

        labelHashSet includePatches;
        includePatches.insert(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 andrayEndFace in global index
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    #include "shootRays.H"

    // Calculate number of visible faces from local index
    labelList nVisibleFaceFaces(nCoarseFaces, 0);

    forAll(rayStartFace, i)
    {
        nVisibleFaceFaces[rayStartFace[i]]++;
    }

    labelListList visibleFaceFaces(nCoarseFaces);

    forAll(nVisibleFaceFaces, facei)
    {
        visibleFaceFaces[facei].setSize(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());

    distributionMap map(globalNumbering, rayEndFace, compactMap);

    labelListIOList IOsubMap
    (
        IOobject
        (
            "subMap",
            mesh.facesInstance(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        map.subMap()
    );
    IOsubMap.write();


    labelListIOList IOconstructMap
    (
        IOobject
        (
            "constructMap",
            mesh.facesInstance(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        map.constructMap()
    );
    IOconstructMap.write();


    IOList<label> consMapDim
    (
        IOobject
        (
            "constructMapDim",
            mesh.facesInstance(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        List<label>(1, map.constructSize())
    );
    consMapDim.write();


    // 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];
        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",
            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, 0.0);

    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().name(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            mesh,
            dimensionedScalar(dimless, 0)
        );

        volScalarField::Boundary& viewFactorFieldBf =
            viewFactorField.boundaryFieldRef();

        label compactI = 0;
        forAll(viewFactorsPatches, i)
        {
            label patchID = viewFactorsPatches[i];
            const labelList& agglom = finalAgglom[patchID];
            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];
                    viewFactorFieldBf[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];

        forAllConstIter(Map<label>, localToCompactMap, iter)
        {
            compactToGlobal[iter()] = globalNumbering.toGlobal
            (
                proci,
                iter.key()
            );
        }
    }


    if (Pstream::master())
    {
        scalarSquareMatrix Fmatrix(totalNCoarseFaces, 0.0);

        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
            ),
            globalFaceFaces
        );
        IOglobalFaceFaces.write();
    }
    else
    {
        labelListList globalFaceFaces(visibleFaceFaces.size());
        forAll(globalFaceFaces, facei)
        {
            globalFaceFaces[facei] = renumber
            (
                compactToGlobal,
                visibleFaceFaces[facei]
            );
        }

        labelListIOList IOglobalFaceFaces
        (
            IOobject
            (
                "globalFaceFaces",
                mesh.facesInstance(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            globalFaceFaces
        );

        IOglobalFaceFaces.write();
    }

    Info<< "End\n" << endl;
    return 0;
}


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