domainDecomposition.C

Go to the documentation of this file.

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

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

#include "domainDecomposition.H"
#include "decompositionMethod.H"
#include "IOobjectList.H"
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
#include "hexRef8Data.H"
#include "cyclicFvPatch.H"
#include "processorCyclicFvPatch.H"
#include "nonConformalCyclicFvPatch.H"
#include "nonConformalProcessorCyclicFvPatch.H"
#include "nonConformalErrorFvPatch.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(domainDecomposition, 0);
}


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::domainDecomposition::decomposePoints()
{
    for (label proci = 0; proci < nProcs(); proci++)
    {
        fvMesh& procMesh = procMeshes_[proci];

        const label pointsCompare =
            compareInstances
            (
                completeMesh().pointsInstance(),
                procMeshes_[proci].pointsInstance()
            );

        if (pointsCompare == -1)
        {
            procMesh.setPoints
            (
                pointField
                (
                    completeMesh().points(),
                    procPointAddressing_[proci]
                )
            );
        }
    }
}


void Foam::domainDecomposition::reconstructPoints()
{
    const label pointsCompare =
        compareInstances
        (
            completeMesh().pointsInstance(),
            procMeshes_[0].pointsInstance()
        );

    if (pointsCompare == 1)
    {
        pointField completePoints(completeMesh().nPoints());

        for (label proci = 0; proci < nProcs(); proci++)
        {
            const fvMesh& procMesh = procMeshes_[proci];

            completePoints.rmap(procMesh.points(), procPointAddressing_[proci]);
        }

        completeMesh_->setPoints(completePoints);
    }
}


bool Foam::domainDecomposition::completeConformal() const
{
    return completeMesh().conformal();
}


bool Foam::domainDecomposition::procsConformal() const
{
    forAll(procMeshes_, proci)
    {
        if (!procMeshes_[proci].conformal())
        {
            return false;
        }
    }

    return true;
}


void Foam::domainDecomposition::unconform()
{
    // ...
    auto checkNonConformalCoupledPatchOrdering = []
    (
        const labelPair& procs,
        const fvPatch& fvp,
        const fvPatch& nbrFvp,
        const labelList& polyFacesPf,
        const labelList& nbrPolyFacesPf
    )
    {
        if (fvp.size() != nbrFvp.size())
        {
            FatalErrorInFunction
                << "Coupled patches " << fvp.name() << " and "
                << nbrFvp.name() << " are not the same size"
                << exit(FatalError);
        }

        if (fvp.size())
        {
            for (label i = 1; i < fvp.size(); ++ i)
            {
                if
                (
                    polyFacesPf[i - 1] > polyFacesPf[i]
                  ? true
                  : polyFacesPf[i - 1] == polyFacesPf[i]
                  ? nbrPolyFacesPf[i - 1] >= nbrPolyFacesPf[i]
                  : false
                )
                {
                    FatalErrorInFunction
                        << "Coupled patches " << fvp.name()
                        << " and " << nbrFvp.name()
                        << " are not in order";

                    if (procs[0] == procs[1])
                    {
                        FatalErrorInFunction
                            << " on process #" << procs[0];
                    }

                    FatalErrorInFunction
                        << exit(FatalError);
                }
            }
        }

        InfoInFunction
            << "Coupled patches " << fvp.name()
            << " and " << nbrFvp.name()
            << " are in order" << endl;
    };

    // ...
    auto checkNonConformalErrorPatchOrdering = []
    (
        const label& proci,
        const fvPatch& fvp,
        const labelList& polyFacesPf
    )
    {
        if (fvp.size())
        {
            for (label i = 1; i < fvp.size(); ++ i)
            {
                if (polyFacesPf[i - 1] > polyFacesPf[i])
                {
                    FatalErrorInFunction
                        << "Error patch " << fvp.name()
                        << " is not in order";

                    if (proci > 0)
                    {
                        FatalErrorInFunction
                            << " on process #" << proci;
                    }

                    FatalErrorInFunction
                        << exit(FatalError);
                }
            }
        }

        InfoInFunction
            << "Error patch " << fvp.name()
            << " is in order" << endl;
    };

    if (completeConformal() && procsConformal())
    {
        // Nothing to do
    }
    else if (!completeConformal() && procsConformal())
    {
        // Construct the reverse of proc-face-face addressing
        labelList faceProcFace(completeMesh().nFaces(), -labelMax);
        forAll(procMeshes_, proci)
        {
            forAll(procFaceAddressing()[proci], procFacei)
            {
                const label facei =
                    mag(procFaceAddressing()[proci][procFacei]) - 1;

                faceProcFace[facei] =
                    faceProcFace[facei] == -labelMax ? procFacei : -1;
            }
        }

        // Map the polyFacesBf from the complete to the processor meshes and
        // unconform the processor meshes. Set dummy data for the face
        // geometry. This should not be used during decomposition.
        forAll(procMeshes_, proci)
        {
            fvMesh& procMesh = procMeshes_[proci];

            const surfaceLabelField::Boundary& faceAddressingBf =
                procFaceAddressingBf()[proci];

            surfaceLabelField::Boundary polyFacesBf
            (
                surfaceLabelField::null(),
                procMesh.polyFacesBf()
            );
            surfaceVectorField Sf(procMesh.Sf().cloneUnSliced());
            surfaceVectorField Cf(procMesh.Cf().cloneUnSliced());

            forAll(procMesh.boundary(), procNccPatchi)
            {
                const fvPatch& fvp = procMesh.boundary()[procNccPatchi];

                if (isA<nonConformalFvPatch>(fvp))
                {
                    const label completeNccPatchi =
                        isA<processorCyclicFvPatch>(fvp)
                      ? refCast<const processorCyclicFvPatch>(fvp)
                       .referPatchID()
                      : procNccPatchi;

                    polyFacesBf[procNccPatchi] =
                        labelField
                        (
                            faceProcFace,
                            labelField
                            (
                                completeMesh().polyFacesBf()[completeNccPatchi],
                                mag(faceAddressingBf[procNccPatchi]) - 1
                            )
                        );

                    const label size = polyFacesBf[procNccPatchi].size();

                    Sf.boundaryFieldRef()[procNccPatchi].resize(size, Zero);
                    Cf.boundaryFieldRef()[procNccPatchi].resize(size, Zero);
                }
            }

            procMesh.unconform
            (
                polyFacesBf,
                Sf,
                Cf,
                NullObjectRef<surfaceScalarField>(),
                false
            );
        }

        // Check ordering
        if (debug)
        {
            forAll(procMeshes_, proci)
            {
                forAll(procMeshes_[proci].boundary(), patchi)
                {
                    const fvPatch& fvp =
                        procMeshes_[proci].boundary()[patchi];

                    // Coupled patches
                    if
                    (
                        isA<nonConformalCoupledFvPatch>(fvp)
                     && refCast<const nonConformalCoupledFvPatch>(fvp).owner()
                    )
                    {
                        const label nccPatchi = patchi;

                        label nbrProci = -1, nbrNccPatchi = -1;
                        if (isA<cyclicFvPatch>(fvp))
                        {
                            nbrProci = proci;
                            nbrNccPatchi =
                                refCast<const cyclicFvPatch>(fvp).nbrPatchID();
                        }
                        else if (isA<processorCyclicFvPatch>(fvp))
                        {
                            typedef processorCyclicFvPatch PcFvp;

                            const PcFvp& pcFvp = refCast<const PcFvp>(fvp);

                            nbrProci = pcFvp.neighbProcNo();

                            const fvBoundaryMesh& nbrFvPatches =
                                procMeshes_[nbrProci].boundary();

                            forAll(nbrFvPatches, nbrNccPatchj)
                            {
                                const fvPatch& nbrFvp =
                                    nbrFvPatches[nbrNccPatchj];

                                if (isA<PcFvp>(nbrFvp))
                                {
                                    const PcFvp& nbrPcFvp =
                                        refCast<const PcFvp>(nbrFvp);

                                    if
                                    (
                                        nbrPcFvp.neighbProcNo()
                                     == proci
                                     && nbrPcFvp.referPatchID()
                                     == pcFvp.referPatch().nbrPatchID()
                                    )
                                    {
                                        nbrNccPatchi = nbrNccPatchj;
                                        break;
                                    }
                                }
                            }

                            if (nbrNccPatchi == -1)
                            {
                                FatalErrorInFunction
                                    << "Opposite processor patch not found for "
                                    << "patch " << fvp.name() << " on proc #"
                                    << proci << exit(FatalError);
                            }
                        }
                        else
                        {
                            FatalErrorInFunction
                                << "Non-conformal-coupled type not recognised "
                                << "for patch " << fvp.name() << " on proc #"
                                << proci << exit(FatalError);
                        }

                        checkNonConformalCoupledPatchOrdering
                        (
                            {proci, nbrProci},
                            procMeshes_[proci].boundary()[nccPatchi],
                            procMeshes_[nbrProci].boundary()[nbrNccPatchi],
                            procMeshes_[proci].polyFacesBf()[nccPatchi],
                            procMeshes_[nbrProci].polyFacesBf()[nbrNccPatchi]
                        );
                    }

                    // Error patches
                    if (isA<nonConformalErrorFvPatch>(fvp))
                    {
                        const label ncePatchi = patchi;

                        checkNonConformalErrorPatchOrdering
                        (
                            proci,
                            procMeshes_[proci].boundary()[ncePatchi],
                            procMeshes_[proci].polyFacesBf()[ncePatchi]
                        );
                    }
                }
            }
        }
    }
    else if (completeConformal() && !procsConformal())
    {
        // Map the polyFacesBf from the processor to the complete meshes and
        // unconform the complete mesh. Set dummy data for the face
        // geometry. This should not be used during decomposition.
        surfaceLabelField::Boundary polyFacesBf
        (
            surfaceLabelField::null(),
            completeMesh().polyFacesBf()
        );
        surfaceVectorField Sf(completeMesh().Sf().cloneUnSliced());
        surfaceVectorField Cf(completeMesh().Cf().cloneUnSliced());

        forAll(procMeshes_, proci)
        {
            const fvMesh& procMesh = procMeshes_[proci];

            const surfaceLabelField::Boundary& faceAddressingBf =
                procFaceAddressingBf()[proci];

            forAll(procMesh.boundary(), procNccPatchi)
            {
                const fvPatch& fvp = procMesh.boundary()[procNccPatchi];

                if (isA<nonConformalFvPatch>(fvp))
                {
                    const label completeNccPatchi =
                        isA<processorCyclicFvPatch>(fvp)
                      ? refCast<const processorCyclicFvPatch>(fvp)
                       .referPatchID()
                      : procNccPatchi;

                    const label size =
                        max
                        (
                            max(mag(faceAddressingBf[procNccPatchi])),
                            polyFacesBf[completeNccPatchi].size()
                        );

                    polyFacesBf[completeNccPatchi].resize(size, -1);
                    Sf.boundaryFieldRef()[completeNccPatchi].resize(size, Zero);
                    Cf.boundaryFieldRef()[completeNccPatchi].resize(size, Zero);

                    polyFacesBf[completeNccPatchi].labelField::rmap
                    (
                        mag
                        (
                            labelField
                            (
                                procFaceAddressing_[proci],
                                procMesh.polyFacesBf()[procNccPatchi]
                            )
                        ) - 1,
                        mag(faceAddressingBf[procNccPatchi]) - 1
                    );
                }
            }
        }

        completeMesh_->unconform
        (
            polyFacesBf,
            Sf,
            Cf,
            NullObjectRef<surfaceScalarField>(),
            false
        );

        // Check ordering
        if (debug)
        {
            forAll(completeMesh().boundary(), patchi)
            {
                const fvPatch& fvp = completeMesh().boundary()[patchi];

                // Coupled patches
                if
                (
                    isA<nonConformalCyclicFvPatch>(fvp)
                 && refCast<const nonConformalCoupledFvPatch>(fvp).owner()
                )
                {
                    const label nccPatchi = patchi;

                    const label nbrNccPatchi =
                        refCast<const nonConformalCyclicFvPatch>(fvp)
                       .nbrPatchID();

                    checkNonConformalCoupledPatchOrdering
                    (
                        {-labelMax, labelMax},
                        completeMesh().boundary()[nccPatchi],
                        completeMesh().boundary()[nbrNccPatchi],
                        completeMesh().polyFacesBf()[nccPatchi],
                        completeMesh().polyFacesBf()[nbrNccPatchi]
                    );
                }

                // Error patches
                if (isA<nonConformalErrorFvPatch>(fvp))
                {
                    const label ncePatchi = patchi;

                    checkNonConformalErrorPatchOrdering
                    (
                        -labelMax,
                        completeMesh().boundary()[ncePatchi],
                        completeMesh().polyFacesBf()[ncePatchi]
                    );
                }
            }
        }
    }
    else // if (!completeConformal() && !procsConformal())
    {
        // Assume everything is consistent and do nothing
    }
}


Foam::label Foam::domainDecomposition::compareInstances
(
    const fileName& a,
    const fileName& b
) const
{
    const word& constant = runTimes_.completeTime().constant();

    if (a == constant && b == constant) return 0;

    if (a == constant) return +1;

    if (b == constant) return -1;

    const scalar aValue = instant(a).value();
    const scalar bValue = instant(b).value();

    if (aValue < bValue) return +1;

    if (aValue > bValue) return -1;

    return 0;
}


void Foam::domainDecomposition::validateComplete() const
{
    if (!completeMesh_.valid())
    {
        FatalErrorInFunction
            << "Complete data requested but complete mesh has not been "
            << "generated or read" << exit(FatalError);
    }
}


void Foam::domainDecomposition::validateProcs() const
{
    if (!procMeshes_.set(0))
    {
        FatalErrorInFunction
            << "Decomposed data requested but decomposed mesh has not been "
            << "generated or read" << exit(FatalError);
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::domainDecomposition::domainDecomposition
(
    const processorRunTimes& runTimes,
    const word& regionName
)
:
    runTimes_(runTimes),
    regionName_(regionName),
    completeMesh_(nullptr),
    procMeshes_(nProcs()),
    procPointAddressing_(nProcs()),
    procFaceAddressing_(nProcs()),
    procCellAddressing_(nProcs()),
    procFaceAddressingBf_()
{}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::domainDecomposition::~domainDecomposition()
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

void Foam::domainDecomposition::readComplete()
{
    completeMesh_.reset
    (
        new fvMesh
        (
            IOobject
            (
                regionName_,
                runTimes_.completeTime().timeName(),
                runTimes_.completeTime(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            false
        )
    );
}


void Foam::domainDecomposition::readProcs()
{
    for (label proci = 0; proci < nProcs(); proci++)
    {
        procMeshes_.set
        (
            proci,
            new fvMesh
            (
                IOobject
                (
                    regionName_,
                    runTimes_.procTimes()[proci].timeName(),
                    runTimes_.procTimes()[proci],
                    IOobject::NO_READ,
                    IOobject::NO_WRITE,
                    false
                ),
                false
            )
        );
    }
}


void Foam::domainDecomposition::readAddressing()
{
    for (label proci = 0; proci < nProcs(); proci++)
    {
        const fvMesh& procMesh = procMeshes_[proci];

        procPointAddressing_[proci] =
            labelIOList
            (
                IOobject
                (
                    "pointProcAddressing",
                    procMesh.facesInstance(),
                    procMesh.meshSubDir,
                    procMesh,
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE,
                    false
                )
            );

        procFaceAddressing_[proci] =
            labelIOList
            (
                IOobject
                (
                    "faceProcAddressing",
                    procMesh.facesInstance(),
                    procMesh.meshSubDir,
                    procMesh,
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE,
                    false
                )
            );

        procCellAddressing_[proci] =
            labelIOList
            (
                IOobject
                (
                    "cellProcAddressing",
                    procMesh.facesInstance(),
                    procMesh.meshSubDir,
                    procMesh,
                    IOobject::MUST_READ,
                    IOobject::NO_WRITE,
                    false
                )
            );
    }
}


Foam::fvMesh::readUpdateState Foam::domainDecomposition::readUpdate()
{
    validateComplete();
    validateProcs();

    // Do read-update on all meshes
    fvMesh::readUpdateState stat =
        completeMesh_->readUpdate(fvMesh::stitchType::nonGeometric);
    forAll(runTimes_.procTimes(), proci)
    {
        fvMesh::readUpdateState procStat =
            procMeshes_[proci].readUpdate(fvMesh::stitchType::nonGeometric);
        if (procStat > stat)
        {
            stat = procStat;
        }
    }

    const label facesCompare =
        compareInstances
        (
            completeMesh().facesInstance(),
            procMeshes_[0].facesInstance()
        );

    const label pointsCompare =
        compareInstances
        (
            completeMesh().pointsInstance(),
            procMeshes_[0].pointsInstance()
        );

    // If the complete mesh has newer topology then we need to decompose
    if (facesCompare == -1)
    {
        decompose();
    }

    // If the processor meshes have newer topology then we need to reconstruct
    if (facesCompare == +1)
    {
        reconstruct();
    }

    // If there has been matching topology change then re-read the addressing
    if (facesCompare == 0 && stat >= fvMesh::TOPO_CHANGE)
    {
        procFaceAddressingBf_.clear();
        readAddressing();
    }

    // If non conformal meshes have changed then clear the finite volume
    // addressing
    if
    (
        (!completeConformal() || !procsConformal())
     && stat != fvMesh::UNCHANGED
    )
    {
        procFaceAddressingBf_.clear();
    }

    // If the complete mesh has newer points then we need to update the
    // processor points and we might need to re-unconform the processor meshes
    if (pointsCompare == -1)
    {
        if (!procsConformal())
        {
            forAll(procMeshes_, proci) procMeshes_[proci].conform();
            decomposePoints();
            unconform();
        }
        else
        {
            decomposePoints();
        }
    }

    // If the processor meshes have newer points then we need to update the
    // complete points and we might need to re-unconform the complete mesh
    if (pointsCompare == +1)
    {
        if (!completeConformal())
        {
            completeMesh_->conform();
            reconstructPoints();
            unconform();
        }
        else
        {
            reconstructPoints();
        }
    }

    return stat;
}


Foam::labelList Foam::domainDecomposition::cellToProc() const
{
    labelList result(completeMesh().nCells());

    forAll(procCellAddressing_, proci)
    {
        forAll(procCellAddressing_[proci], procCelli)
        {
            result[procCellAddressing_[proci][procCelli]] = proci;
        }
    }

    return result;
}


const Foam::PtrList<Foam::surfaceLabelField::Boundary>&
Foam::domainDecomposition::procFaceAddressingBf() const
{
    validateComplete();
    validateProcs();

    if (procFaceAddressingBf_.empty())
    {
        // Construct cell-to-proc addressing (this is a bit wasteful, as it been
        // available before but was thrown away)
        const labelList cellToProc(this->cellToProc());

        // Map from reference patch and processors to the interface patch
        typedef HashTable<label, labelPair, Hash<labelPair>> labelPairTable;
        List<labelPairTable> refPatchProcPatchTable
        (
            completeMesh().boundary().size()
        );
        forAll(procMeshes_, proci)
        {
            const fvMesh& procMesh = procMeshes_[proci];

            forAll(procMesh.boundary(), procPatchi)
            {
                const fvPatch& fvp = procMesh.boundary()[procPatchi];

                if (isA<cyclicFvPatch>(fvp))
                {
                    refPatchProcPatchTable[procPatchi].insert
                    (
                        labelPair(proci, proci),
                        procPatchi
                    );
                }
                else if (isA<processorCyclicFvPatch>(fvp))
                {
                    const processorCyclicFvPatch& pcFvp =
                        refCast<const processorCyclicFvPatch>(fvp);

                    refPatchProcPatchTable[pcFvp.referPatchID()].insert
                    (
                        labelPair(proci, pcFvp.neighbProcNo()),
                        procPatchi
                    );
                }
            }
        }

        // Build non-conformal finite volume face addressing for each processor
        List<List<DynamicList<label>>>
            nonConformalProcFaceAddressingBf(nProcs());
        forAll(nonConformalProcFaceAddressingBf, proci)
        {
            nonConformalProcFaceAddressingBf[proci].resize
            (
                procMeshes_[proci].boundary().size()
            );
        }
        if (completeMesh().conformal() && procMeshes_[0].conformal())
        {
            // Nothing to do
        }
        else if (!completeMesh().conformal())
        {
            // Decompose non-conformal addressing
            const surfaceLabelField::Boundary& polyFacesBf =
                completeMesh().polyFacesBf();

            // Cyclic patches
            forAll(completeMesh().boundary(), nccPatchi)
            {
                const fvPatch& fvp = completeMesh().boundary()[nccPatchi];

                if (!isA<nonConformalCyclicFvPatch>(fvp)) continue;

                const nonConformalCyclicFvPatch& nccFvp =
                    refCast<const nonConformalCyclicFvPatch>(fvp);

                if (!nccFvp.owner()) continue;

                const label nccNbrPatchi = nccFvp.nbrPatchID();

                forAll(polyFacesBf[nccPatchi], nccPatchFacei)
                {
                    const label facei = polyFacesBf[nccPatchi][nccPatchFacei];
                    const label celli = completeMesh().faceOwner()[facei];
                    const label proci = cellToProc[celli];

                    const label nbrFacei =
                        polyFacesBf[nccNbrPatchi][nccPatchFacei];
                    const label nbrCelli =
                        completeMesh().faceOwner()[nbrFacei];
                    const label nbrProci = cellToProc[nbrCelli];

                    const label procNccPatchi =
                        refPatchProcPatchTable
                        [nccPatchi][labelPair(proci, nbrProci)];
                    const label nbrProcNccPatchi =
                        refPatchProcPatchTable
                        [nccNbrPatchi][labelPair(nbrProci, proci)];

                    nonConformalProcFaceAddressingBf[proci][procNccPatchi]
                       .append(nccPatchFacei + 1);
                    nonConformalProcFaceAddressingBf[nbrProci][nbrProcNccPatchi]
                       .append(nccPatchFacei + 1);
                }
            }

            // Error patches
            forAll(completeMesh().boundary(), ncePatchi)
            {
                const fvPatch& fvp = completeMesh().boundary()[ncePatchi];

                if (!isA<nonConformalErrorFvPatch>(fvp)) continue;

                forAll(polyFacesBf[ncePatchi], ncePatchFacei)
                {
                    const label facei = polyFacesBf[ncePatchi][ncePatchFacei];
                    const label celli = completeMesh().faceOwner()[facei];
                    const label proci = cellToProc[celli];

                    nonConformalProcFaceAddressingBf[proci][ncePatchi]
                       .append(ncePatchFacei + 1);
                }
            }
        }
        else // if (!procMeshes_[0].conformal())
        {
            // Reconstruct non-conformal addressing

            // Cyclic patches
            forAll(completeMesh().boundary(), nccPatchi)
            {
                const fvPatch& fvp = completeMesh().boundary()[nccPatchi];

                if (!isA<nonConformalCyclicFvPatch>(fvp)) continue;

                const nonConformalCyclicFvPatch& nccFvp =
                    refCast<const nonConformalCyclicFvPatch>(fvp);

                if (!nccFvp.owner()) continue;

                const label nccNbrPatchi = nccFvp.nbrPatchID();

                label count = 0;
                labelPairTable procCounts;
                forAllConstIter
                (
                    labelPairTable,
                    refPatchProcPatchTable[nccPatchi],
                    iter
                )
                {
                    procCounts.insert(iter.key(), 0);
                }

                while (true)
                {
                    labelPair procNbrProc(labelMax, labelMax);
                    labelPair faceNbrFace(labelMax, labelMax);

                    forAllConstIter(labelPairTable, procCounts, iter)
                    {
                        const label proci = iter.key().first();
                        const label nbrProci = iter.key().second();

                        const labelPair procNbrProcStar(proci, nbrProci);
                        const labelPair nbrProcProcStar(nbrProci, proci);

                        const label procNccPatchi =
                            refPatchProcPatchTable
                            [nccPatchi][procNbrProcStar];
                        const label nbrProcNccPatchi =
                            refPatchProcPatchTable
                            [nccNbrPatchi][nbrProcProcStar];

                        const label size =
                            procMeshes_[proci]
                           .polyFacesBf()[procNccPatchi]
                           .size();

                        if (iter() >= size) continue;

                        const label procFacei =
                            procMeshes_[proci].polyFacesBf()
                            [procNccPatchi][iter()];
                        const label nbrProcFacei =
                            procMeshes_[nbrProci].polyFacesBf()
                            [nbrProcNccPatchi][iter()];

                        const labelPair faceNbrFaceStar
                        (
                            procFaceAddressing_[proci][procFacei] - 1,
                            procFaceAddressing_[nbrProci][nbrProcFacei] - 1
                        );

                        if (faceNbrFace > faceNbrFaceStar)
                        {
                            procNbrProc = procNbrProcStar;
                            faceNbrFace = faceNbrFaceStar;
                        }
                    }

                    if (faceNbrFace == labelPair(labelMax, labelMax))
                    {
                        break;
                    }
                    else
                    {
                        const label proci = procNbrProc.first();
                        const label nbrProci = procNbrProc.second();

                        const labelPair nbrProcProc(nbrProci, proci);

                        const label procNccPatchi =
                            refPatchProcPatchTable[nccPatchi][procNbrProc];
                        const label nbrProcNccPatchi =
                            refPatchProcPatchTable[nccNbrPatchi][nbrProcProc];

                        nonConformalProcFaceAddressingBf
                            [proci][procNccPatchi].append(count + 1);
                        nonConformalProcFaceAddressingBf
                            [nbrProci][nbrProcNccPatchi].append(count + 1);

                        count ++;
                        procCounts[procNbrProc] ++;
                    }
                }
            }

            // Error patches
            forAll(completeMesh().boundary(), ncePatchi)
            {
                const fvPatch& fvp = completeMesh().boundary()[ncePatchi];

                if (!isA<nonConformalErrorFvPatch>(fvp)) continue;

                label count = 0;
                labelList procCounts(nProcs(), 0);

                while (true)
                {
                    label facei = labelMax, proci = labelMax;

                    forAll(procCounts, procStari)
                    {
                        const label size =
                            procMeshes_[procStari]
                           .polyFacesBf()[ncePatchi]
                           .size();

                        if (procCounts[procStari] >= size) continue;

                        const label procFacei =
                            procMeshes_[procStari].polyFacesBf()
                            [ncePatchi][procCounts[procStari]];

                        const label faceStari =
                            procFaceAddressing_[procStari][procFacei] - 1;

                        if (facei > faceStari)
                        {
                            facei = faceStari;
                            proci = procStari;
                        }
                    }

                    if (facei == labelMax)
                    {
                        break;
                    }
                    else
                    {
                        nonConformalProcFaceAddressingBf
                            [proci][ncePatchi].append(count + 1);

                        count ++;
                        procCounts[proci] ++;
                    }
                }
            }
        }

        // Build finite volume face addressing boundary fields
        procFaceAddressingBf_.resize(nProcs());
        forAll(procMeshes_, proci)
        {
            const fvMesh& procMesh = procMeshes_[proci];

            procFaceAddressingBf_.set
            (
                proci,
                new surfaceLabelField::Boundary
                (
                    procMesh.boundary(),
                    surfaceLabelField::null(),
                    calculatedFvsPatchLabelField::typeName
                )
            );

            forAll(procMesh.boundary(), procPatchi)
            {
                const fvPatch& fvp = procMesh.boundary()[procPatchi];

                if (isA<nonConformalFvPatch>(fvp))
                {
                    procFaceAddressingBf_[proci][procPatchi] =
                        nonConformalProcFaceAddressingBf[proci][procPatchi];
                }
                else if (isA<processorCyclicFvPatch>(fvp))
                {
                    const label completePatchi =
                        refCast<const processorCyclicFvPatch>(fvp)
                       .referPatchID();

                    procFaceAddressingBf_[proci][procPatchi] =
                        mag(fvp.patchSlice(procFaceAddressing_[proci]))
                      - completeMesh().boundaryMesh()[completePatchi].start();
                }
                else if (isA<processorFvPatch>(fvp))
                {
                    procFaceAddressingBf_[proci][procPatchi] =
                        fvp.patchSlice(procFaceAddressing_[proci]);
                }
                else
                {
                    procFaceAddressingBf_[proci][procPatchi] =
                        mag(fvp.patchSlice(procFaceAddressing_[proci]))
                      - completeMesh().boundaryMesh()[procPatchi].start();
                }
            }
        }
    }

    return procFaceAddressingBf_;
}


void Foam::domainDecomposition::writeAddressing() const
{
    for (label proci = 0; proci < nProcs(); proci++)
    {
        const fvMesh& procMesh = procMeshes_[proci];

        labelIOList pointProcAddressing
        (
            IOobject
            (
                "pointProcAddressing",
                procMesh.facesInstance(),
                procMesh.meshSubDir,
                procMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            procPointAddressing_[proci]
        );
        pointProcAddressing.write();

        labelIOList faceProcAddressing
        (
            IOobject
            (
                "faceProcAddressing",
                procMesh.facesInstance(),
                procMesh.meshSubDir,
                procMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            procFaceAddressing_[proci]
        );
        faceProcAddressing.write();

        labelIOList cellProcAddressing
        (
            IOobject
            (
                "cellProcAddressing",
                procMesh.facesInstance(),
                procMesh.meshSubDir,
                procMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            procCellAddressing_[proci]
        );
        cellProcAddressing.write();
    }
}


void Foam::domainDecomposition::writeComplete(const bool doSets) const
{
    // Set the precision of the points data to be min 10
    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));

    // Write the processor mesh
    completeMesh().write();

    // Sets
    if (doSets)
    {
        // Scan to find all sets
        HashTable<label> cSetNames;
        HashTable<label> fSetNames;
        HashTable<label> pSetNames;
        for (label proci = 0; proci < nProcs(); proci++)
        {
            const fvMesh& procMesh = procMeshes_[proci];

            IOobjectList objects
            (
                procMesh,
                runTimes_.procTimes()[proci].timeName(),
                polyMesh::meshSubDir/"sets"
            );

            IOobjectList cSets(objects.lookupClass(cellSet::typeName));
            forAllConstIter(IOobjectList, cSets, iter)
            {
                cSetNames.insert(iter.key(), cSetNames.size());
            }

            IOobjectList fSets(objects.lookupClass(faceSet::typeName));
            forAllConstIter(IOobjectList, fSets, iter)
            {
                fSetNames.insert(iter.key(), fSetNames.size());
            }
            IOobjectList pSets(objects.lookupClass(pointSet::typeName));
            forAllConstIter(IOobjectList, pSets, iter)
            {
                pSetNames.insert(iter.key(), pSetNames.size());
            }
        }

        // Reconstruct all sets
        if (cSetNames.size() || fSetNames.size() || pSetNames.size())
        {
            PtrList<cellSet> cellSets(cSetNames.size());
            PtrList<faceSet> faceSets(fSetNames.size());
            PtrList<pointSet> pointSets(pSetNames.size());

            Info<< "Reconstructing sets:" << endl;
            if (cSetNames.size())
            {
                Info<< "    cellSets "
                    << cSetNames.sortedToc() << endl;
            }
            if (fSetNames.size())
            {
                Info<< "    faceSets "
                    << fSetNames.sortedToc() << endl;
            }
            if (pSetNames.size())
            {
                Info<< "    pointSets "
                    << pSetNames.sortedToc() << endl;
            }

            // Load sets
            for (label proci = 0; proci < nProcs(); proci++)
            {
                const fvMesh& procMesh = procMeshes_[proci];

                IOobjectList objects
                (
                    procMesh,
                    runTimes_.procTimes()[proci].timeName(),
                    polyMesh::meshSubDir/"sets"
                );

                // cellSets
                const labelList& cellMap = procCellAddressing()[proci];

                IOobjectList cSets
                (
                    objects.lookupClass(cellSet::typeName)
                );

                forAllConstIter(IOobjectList, cSets, iter)
                {
                    // Load cellSet
                    const cellSet procSet(*iter());
                    label setI = cSetNames[iter.key()];
                    if (!cellSets.set(setI))
                    {
                        cellSets.set
                        (
                            setI,
                            new cellSet
                            (
                                completeMesh(),
                                iter.key(),
                                procSet.size()
                            )
                        );
                    }
                    cellSet& cSet = cellSets[setI];
                    cSet.instance() = runTimes_.completeTime().timeName();

                    forAllConstIter(cellSet, procSet, iter)
                    {
                        cSet.insert(cellMap[iter.key()]);
                    }
                }

                // faceSets
                const labelList& faceMap = procFaceAddressing()[proci];

                IOobjectList fSets
                (
                    objects.lookupClass(faceSet::typeName)
                );

                forAllConstIter(IOobjectList, fSets, iter)
                {
                    // Load faceSet
                    const faceSet procSet(*iter());
                    label setI = fSetNames[iter.key()];
                    if (!faceSets.set(setI))
                    {
                        faceSets.set
                        (
                            setI,
                            new faceSet
                            (
                                completeMesh(),
                                iter.key(),
                                procSet.size()
                            )
                        );
                    }
                    faceSet& fSet = faceSets[setI];
                    fSet.instance() = runTimes_.completeTime().timeName();

                    forAllConstIter(faceSet, procSet, iter)
                    {
                        fSet.insert(mag(faceMap[iter.key()]) - 1);
                    }
                }

                // pointSets
                const labelList& pointMap = procPointAddressing()[proci];

                IOobjectList pSets
                (
                    objects.lookupClass(pointSet::typeName)
                );

                forAllConstIter(IOobjectList, pSets, iter)
                {
                    // Load pointSet
                    const pointSet propSet(*iter());
                    label setI = pSetNames[iter.key()];
                    if (!pointSets.set(setI))
                    {
                        pointSets.set
                        (
                            setI,
                            new pointSet
                            (
                                completeMesh(),
                                iter.key(),
                                propSet.size()
                            )
                        );
                    }
                    pointSet& pSet = pointSets[setI];
                    pSet.instance() = runTimes_.completeTime().timeName();

                    forAllConstIter(pointSet, propSet, iter)
                    {
                        pSet.insert(pointMap[iter.key()]);
                    }
                }
            }

            // Write sets
            forAll(cellSets, i)
            {
                cellSets[i].write();
            }
            forAll(faceSets, i)
            {
                faceSets[i].write();
            }
            forAll(pointSets, i)
            {
                pointSets[i].write();
            }
        }
    }

    // Refinement data (if any)
    UPtrList<const labelList> cellMaps(nProcs());
    UPtrList<const labelList> pointMaps(nProcs());
    PtrList<const hexRef8Data> refinementDatas(nProcs());
    for (label proci = 0; proci < nProcs(); proci++)
    {
        const fvMesh& procMesh = procMeshes_[proci];

        cellMaps.set(proci, &procCellAddressing_[proci]);
        pointMaps.set(proci, &procPointAddressing_[proci]);
        refinementDatas.set
        (
            proci,
            new hexRef8Data
            (
                IOobject
                (
                    "dummy",
                    completeMesh().facesInstance(),
                    polyMesh::meshSubDir,
                    procMesh,
                    IOobject::READ_IF_PRESENT,
                    IOobject::NO_WRITE,
                    false
                )
            )
        );
    }
    hexRef8Data
    (
        IOobject
        (
            "dummy",
            completeMesh().facesInstance(),
            polyMesh::meshSubDir,
            completeMesh(),
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        cellMaps,
        pointMaps,
        refinementDatas
    ).write();

    // Write decomposition addressing
    writeAddressing();
}


void Foam::domainDecomposition::writeProcs(const bool doSets) const
{
    // Read sets
    PtrList<const cellSet> cellSets;
    PtrList<const faceSet> faceSets;
    PtrList<const pointSet> pointSets;
    if (doSets)
    {
        IOobjectList objects
        (
            completeMesh(),
            completeMesh().facesInstance(),
            "polyMesh/sets"
        );
        {
            IOobjectList cSets(objects.lookupClass(cellSet::typeName));
            forAllConstIter(IOobjectList, cSets, iter)
            {
                cellSets.append(new cellSet(*iter()));
            }
        }
        {
            IOobjectList fSets(objects.lookupClass(faceSet::typeName));
            forAllConstIter(IOobjectList, fSets, iter)
            {
                faceSets.append(new faceSet(*iter()));
            }
        }
        {
            IOobjectList pSets(objects.lookupClass(pointSet::typeName));
            forAllConstIter(IOobjectList, pSets, iter)
            {
                pointSets.append(new pointSet(*iter()));
            }
        }
    }

    // Read refinement data (if any)
    hexRef8Data refinementData
    (
        IOobject
        (
            "dummy",
            completeMesh().facesInstance(),
            polyMesh::meshSubDir,
            completeMesh(),
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        )
    );

    // Write out the meshes
    for (label proci = 0; proci < nProcs(); proci++)
    {
        const fvMesh& procMesh = procMeshes_[proci];

        // Set the precision of the points data to be min 10
        IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));

        // Write the processor mesh
        procMesh.write();

        // Write any sets
        if (doSets)
        {
            forAll(cellSets, i)
            {
                const cellSet& cs = cellSets[i];
                cellSet set(procMesh, cs.name(), cs.size()/nProcs());
                forAll(procCellAddressing_[proci], i)
                {
                    if (cs.found(procCellAddressing_[proci][i]))
                    {
                        set.insert(i);
                    }
                }
                set.write();
            }
            forAll(faceSets, i)
            {
                const faceSet& cs = faceSets[i];
                faceSet set(procMesh, cs.name(), cs.size()/nProcs());
                forAll(procFaceAddressing_[proci], i)
                {
                    if (cs.found(mag(procFaceAddressing_[proci][i])-1))
                    {
                        set.insert(i);
                    }
                }
                set.write();
            }
            forAll(pointSets, i)
            {
                const pointSet& cs = pointSets[i];
                pointSet set(procMesh, cs.name(), cs.size()/nProcs());
                forAll(procPointAddressing_[proci], i)
                {
                    if (cs.found(procPointAddressing_[proci][i]))
                    {
                        set.insert(i);
                    }
                }
                set.write();
            }
        }

        // Write refinement data (if any)
        hexRef8Data
        (
            IOobject
            (
                "dummy",
                completeMesh().facesInstance(),
                polyMesh::meshSubDir,
                procMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            refinementData,
            procCellAddressing_[proci],
            procPointAddressing_[proci]
        ).write();
    }

    // Write decomposition addressing
    writeAddressing();
}


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