domainDecomposition.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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 "dictionary.H"
#include "labelIOList.H"
#include "processorPolyPatch.H"
#include "processorCyclicPolyPatch.H"
#include "fvMesh.H"
#include "OSspecific.H"
#include "Map.H"
#include "DynamicList.H"
#include "fvFieldDecomposer.H"
#include "IOobjectList.H"
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
#include "decompositionModel.H"
#include "hexRef8Data.H"
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
void Foam::domainDecomposition::mark
(
    const labelList& zoneElems,
    const label zoneI,
    labelList& elementToZone
)
{
    forAll(zoneElems, i)
    {
        label pointi = zoneElems[i];
 
        if (elementToZone[pointi] == -1)
        {
            // First occurrence
            elementToZone[pointi] = zoneI;
        }
        else if (elementToZone[pointi] >= 0)
        {
            // Multiple zones
            elementToZone[pointi] = -2;
        }
    }
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
// from components
Foam::domainDecomposition::domainDecomposition
(
    const IOobject& io,
    const fileName& decompDictFile
)
:
    fvMesh(io),
    facesInstancePointsPtr_
    (
        pointsInstance() != facesInstance()
      ? new pointIOField
        (
            IOobject
            (
                "points",
                facesInstance(),
                polyMesh::meshSubDir,
                *this,
                IOobject::MUST_READ,
                IOobject::NO_WRITE,
                false
            )
        )
      : NULL
    ),
    decompDictFile_(decompDictFile),
    nProcs_
    (
        readInt
        (
            decompositionModel::New
            (
                *this,
                decompDictFile
            ).lookup("numberOfSubdomains")
        )
    ),
    distributed_(false),
    cellToProc_(nCells()),
    procPointAddressing_(nProcs_),
    procFaceAddressing_(nProcs_),
    procCellAddressing_(nProcs_),
    procPatchSize_(nProcs_),
    procPatchStartIndex_(nProcs_),
    procNeighbourProcessors_(nProcs_),
    procProcessorPatchSize_(nProcs_),
    procProcessorPatchStartIndex_(nProcs_),
    procProcessorPatchSubPatchIDs_(nProcs_),
    procProcessorPatchSubPatchStarts_(nProcs_)
{
    decompositionModel::New
    (
        *this,
        decompDictFile
    ).readIfPresent("distributed", distributed_);
}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
Foam::domainDecomposition::~domainDecomposition()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
bool Foam::domainDecomposition::writeDecomposition(const bool decomposeSets)
{
    Info<< "\nConstructing processor meshes" << endl;
 
    // Mark point/faces/cells that are in zones.
    // -1   : not in zone
    // -2   : in multiple zones
    // >= 0 : in single given zone
    // This will give direct lookup of elements that are in a single zone
    // and we'll only have to revert back to searching through all zones
    // for the duplicate elements
 
    // Point zones
    labelList pointToZone(points().size(), -1);
 
    forAll(pointZones(), zoneI)
    {
        mark(pointZones()[zoneI], zoneI, pointToZone);
    }
 
    // Face zones
    labelList faceToZone(faces().size(), -1);
 
    forAll(faceZones(), zoneI)
    {
        mark(faceZones()[zoneI], zoneI, faceToZone);
    }
 
    // Cell zones
    labelList cellToZone(nCells(), -1);
 
    forAll(cellZones(), zoneI)
    {
        mark(cellZones()[zoneI], zoneI, cellToZone);
    }
 
 
    PtrList<const cellSet> cellSets;
    PtrList<const faceSet> faceSets;
    PtrList<const pointSet> pointSets;
    if (decomposeSets)
    {
        // Read sets
        IOobjectList objects(*this, 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()));
            }
        }
    }
 
 
    // Load refinement data (if any)
    hexRef8Data baseMeshData
    (
        IOobject
        (
            "dummy",
            facesInstance(),
            polyMesh::meshSubDir,
            *this,
            IOobject::READ_IF_PRESENT,
            IOobject::NO_WRITE,
            false
        )
    );
 
 
 
    label maxProcCells = 0;
    label totProcFaces = 0;
    label maxProcPatches = 0;
    label totProcPatches = 0;
    label maxProcFaces = 0;
 
 
    // Write out the meshes
    for (label procI = 0; procI < nProcs_; procI++)
    {
        // Create processor points
        const labelList& curPointLabels = procPointAddressing_[procI];
 
        const pointField& meshPoints = points();
 
        labelList pointLookup(nPoints(), -1);
 
        pointField procPoints(curPointLabels.size());
 
        forAll(curPointLabels, pointi)
        {
            procPoints[pointi] = meshPoints[curPointLabels[pointi]];
 
            pointLookup[curPointLabels[pointi]] = pointi;
        }
 
        // Create processor faces
        const labelList& curFaceLabels = procFaceAddressing_[procI];
 
        const faceList& meshFaces = faces();
 
        labelList faceLookup(nFaces(), -1);
 
        faceList procFaces(curFaceLabels.size());
 
        forAll(curFaceLabels, facei)
        {
            // Mark the original face as used
            // Remember to decrement the index by one (turning index)
            //
            label curF = mag(curFaceLabels[facei]) - 1;
 
            faceLookup[curF] = facei;
 
            // get the original face
            labelList origFaceLabels;
 
            if (curFaceLabels[facei] >= 0)
            {
                // face not turned
                origFaceLabels = meshFaces[curF];
            }
            else
            {
                origFaceLabels = meshFaces[curF].reverseFace();
            }
 
            // translate face labels into local point list
            face& procFaceLabels = procFaces[facei];
 
            procFaceLabels.setSize(origFaceLabels.size());
 
            forAll(origFaceLabels, pointi)
            {
                procFaceLabels[pointi] = pointLookup[origFaceLabels[pointi]];
            }
        }
 
        // Create processor cells
        const labelList& curCellLabels = procCellAddressing_[procI];
 
        const cellList& meshCells = cells();
 
        cellList procCells(curCellLabels.size());
 
        forAll(curCellLabels, celli)
        {
            const labelList& origCellLabels = meshCells[curCellLabels[celli]];
 
            cell& curCell = procCells[celli];
 
            curCell.setSize(origCellLabels.size());
 
            forAll(origCellLabels, cellFaceI)
            {
                curCell[cellFaceI] = faceLookup[origCellLabels[cellFaceI]];
            }
        }
 
        // Create processor mesh without a boundary
 
        fileName processorCasePath
        (
            time().caseName()/fileName(word("processor") + Foam::name(procI))
        );
 
        // make the processor directory
        mkDir(time().rootPath()/processorCasePath);
 
        // create a database
        Time processorDb
        (
            Time::controlDictName,
            time().rootPath(),
            processorCasePath,
            word("system"),
            word("constant")
        );
        processorDb.setTime(time());
 
        // create the mesh. Two situations:
        // - points and faces come from the same time ('instance'). The mesh
        //   will get constructed in the same instance.
        // - points come from a different time (moving mesh cases).
        //   It will read the points belonging to the faces instance and
        //   construct the procMesh with it which then gets handled as above.
        //   (so with 'old' geometry).
        //   Only at writing time will it additionally write the current
        //   points.
 
        autoPtr<polyMesh> procMeshPtr;
 
        if (facesInstancePointsPtr_.valid())
        {
            // Construct mesh from facesInstance.
            pointField facesInstancePoints
            (
                facesInstancePointsPtr_(),
                curPointLabels
            );
 
            procMeshPtr.reset
            (
                new polyMesh
                (
                    IOobject
                    (
                        this->polyMesh::name(), // region of undecomposed mesh
                        facesInstance(),
                        processorDb
                    ),
                    xferMove(facesInstancePoints),
                    xferMove(procFaces),
                    xferMove(procCells)
                )
            );
        }
        else
        {
            procMeshPtr.reset
            (
                new polyMesh
                (
                    IOobject
                    (
                        this->polyMesh::name(), // region of undecomposed mesh
                        facesInstance(),
                        processorDb
                    ),
                    xferMove(procPoints),
                    xferMove(procFaces),
                    xferMove(procCells)
                )
            );
        }
        polyMesh& procMesh = procMeshPtr();
 
 
        // Create processor boundary patches
        const labelList& curPatchSizes = procPatchSize_[procI];
 
        const labelList& curPatchStarts = procPatchStartIndex_[procI];
 
        const labelList& curNeighbourProcessors =
            procNeighbourProcessors_[procI];
 
        const labelList& curProcessorPatchSizes =
            procProcessorPatchSize_[procI];
 
        const labelList& curProcessorPatchStarts =
            procProcessorPatchStartIndex_[procI];
 
        const labelListList& curSubPatchIDs =
            procProcessorPatchSubPatchIDs_[procI];
 
        const labelListList& curSubStarts =
            procProcessorPatchSubPatchStarts_[procI];
 
        const polyPatchList& meshPatches = boundaryMesh();
 
 
        // Count the number of inter-proc patches
        label nInterProcPatches = 0;
        forAll(curSubPatchIDs, procPatchI)
        {
            //Info<< "For processor " << procI
            //    << " have to destination processor "
            //    << curNeighbourProcessors[procPatchI] << endl;
            //
            //forAll(curSubPatchIDs[procPatchI], i)
            //{
            //    Info<< "    from patch:" << curSubPatchIDs[procPatchI][i]
            //        << " starting at:" << curSubStarts[procPatchI][i]
            //        << endl;
            //}
 
            nInterProcPatches += curSubPatchIDs[procPatchI].size();
        }
 
        //Info<< "For processor " << procI
        //    << " have " << nInterProcPatches
        //    << " patches to neighbouring processors" << endl;
 
 
        List<polyPatch*> procPatches
        (
            curPatchSizes.size()
          + nInterProcPatches,          //curProcessorPatchSizes.size(),
            reinterpret_cast<polyPatch*>(0)
        );
 
        label nPatches = 0;
 
        forAll(curPatchSizes, patchi)
        {
            // Get the face labels consistent with the field mapping
            // (reuse the patch field mappers)
            const polyPatch& meshPatch = meshPatches[patchi];
 
            fvFieldDecomposer::patchFieldDecomposer patchMapper
            (
                SubList<label>
                (
                    curFaceLabels,
                    curPatchSizes[patchi],
                    curPatchStarts[patchi]
                ),
                meshPatch.start()
            );
 
            // Map existing patches
            procPatches[nPatches] = meshPatch.clone
            (
                procMesh.boundaryMesh(),
                nPatches,
                patchMapper.directAddressing(),
                curPatchStarts[patchi]
            ).ptr();
 
            nPatches++;
        }
 
        forAll(curProcessorPatchSizes, procPatchI)
        {
            const labelList& subPatchID = curSubPatchIDs[procPatchI];
            const labelList& subStarts = curSubStarts[procPatchI];
 
            label curStart = curProcessorPatchStarts[procPatchI];
 
            forAll(subPatchID, i)
            {
                label size =
                (
                    i < subPatchID.size()-1
                  ? subStarts[i+1] - subStarts[i]
                  : curProcessorPatchSizes[procPatchI] - subStarts[i]
                );
 
//                Info<< "From processor:" << procI << endl
//                    << "  to processor:" << curNeighbourProcessors[procPatchI]
//                    << endl
//                    << "    via patch:" << subPatchID[i] << endl
//                    << "    start    :" << curStart << endl
//                    << "    size     :" << size << endl;
 
                if (subPatchID[i] == -1)
                {
                    // From internal faces
                    procPatches[nPatches] =
                        new processorPolyPatch
                        (
                            word("procBoundary") + Foam::name(procI)
                          + "to"
                          + Foam::name(curNeighbourProcessors[procPatchI]),
                            size,
                            curStart,
                            nPatches,
                            procMesh.boundaryMesh(),
                            procI,
                            curNeighbourProcessors[procPatchI]
                        );
                }
                else
                {
                    const coupledPolyPatch& pcPatch
                        = refCast<const coupledPolyPatch>
                          (
                              boundaryMesh()[subPatchID[i]]
                          );
 
                    // From cyclic
                    const word& referPatch = pcPatch.name();
 
                    procPatches[nPatches] =
                        new processorCyclicPolyPatch
                        (
                            word("procBoundary") + Foam::name(procI)
                          + "to"
                          + Foam::name(curNeighbourProcessors[procPatchI])
                          + "through"
                          + referPatch,
                            size,
                            curStart,
                            nPatches,
                            procMesh.boundaryMesh(),
                            procI,
                            curNeighbourProcessors[procPatchI],
                            referPatch,
                            pcPatch.transform()
                        );
                }
 
                curStart += size;
 
                nPatches++;
            }
        }
 
 
        //forAll(procPatches, patchI)
        //{
        //    Pout<< "    " << patchI
        //        << '\t' << "name:" << procPatches[patchI]->name()
        //        << '\t' << "type:" << procPatches[patchI]->type()
        //        << '\t' << "size:" << procPatches[patchI]->size()
        //        << endl;
        //}
 
        // Add boundary patches
        procMesh.addPatches(procPatches);
 
        // Create and add zones
 
        // Point zones
        {
            const pointZoneMesh& pz = pointZones();
 
            // Go through all the zoned points and find out if they
            // belong to a zone.  If so, add it to the zone as
            // necessary
            List<DynamicList<label> > zonePoints(pz.size());
 
            // Estimate size
            forAll(zonePoints, zoneI)
            {
                zonePoints[zoneI].setCapacity(pz[zoneI].size() / nProcs_);
            }
 
            // Use the pointToZone map to find out the single zone (if any),
            // use slow search only for shared points.
            forAll(curPointLabels, pointi)
            {
                label curPoint = curPointLabels[pointi];
 
                label zoneI = pointToZone[curPoint];
 
                if (zoneI >= 0)
                {
                    // Single zone.
                    zonePoints[zoneI].append(pointi);
                }
                else if (zoneI == -2)
                {
                    // Multiple zones. Lookup.
                    forAll(pz, zoneI)
                    {
                        label index = pz[zoneI].whichPoint(curPoint);
 
                        if (index != -1)
                        {
                            zonePoints[zoneI].append(pointi);
                        }
                    }
                }
            }
 
            procMesh.pointZones().clearAddressing();
            procMesh.pointZones().setSize(zonePoints.size());
            forAll(zonePoints, zoneI)
            {
                procMesh.pointZones().set
                (
                    zoneI,
                    pz[zoneI].clone
                    (
                        procMesh.pointZones(),
                        zoneI,
                        zonePoints[zoneI].shrink()
                    )
                );
            }
 
            if (pz.size())
            {
                // Force writing on all processors
                procMesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
            }
        }
 
        // Face zones
        {
            const faceZoneMesh& fz = faceZones();
 
            // Go through all the zoned face and find out if they
            // belong to a zone.  If so, add it to the zone as
            // necessary
            List<DynamicList<label> > zoneFaces(fz.size());
            List<DynamicList<bool> > zoneFaceFlips(fz.size());
 
            // Estimate size
            forAll(zoneFaces, zoneI)
            {
                label procSize = fz[zoneI].size() / nProcs_;
 
                zoneFaces[zoneI].setCapacity(procSize);
                zoneFaceFlips[zoneI].setCapacity(procSize);
            }
 
            // Go through all the zoned faces and find out if they
            // belong to a zone.  If so, add it to the zone as
            // necessary
            forAll(curFaceLabels, facei)
            {
                // Remember to decrement the index by one (turning index)
                //
                label curF = mag(curFaceLabels[facei]) - 1;
 
                label zoneI = faceToZone[curF];
 
                if (zoneI >= 0)
                {
                    // Single zone. Add the face
                    zoneFaces[zoneI].append(facei);
 
                    label index = fz[zoneI].whichFace(curF);
 
                    bool flip = fz[zoneI].flipMap()[index];
 
                    if (curFaceLabels[facei] < 0)
                    {
                        flip = !flip;
                    }
 
                    zoneFaceFlips[zoneI].append(flip);
                }
                else if (zoneI == -2)
                {
                    // Multiple zones. Lookup.
                    forAll(fz, zoneI)
                    {
                        label index = fz[zoneI].whichFace(curF);
 
                        if (index != -1)
                        {
                            zoneFaces[zoneI].append(facei);
 
                            bool flip = fz[zoneI].flipMap()[index];
 
                            if (curFaceLabels[facei] < 0)
                            {
                                flip = !flip;
                            }
 
                            zoneFaceFlips[zoneI].append(flip);
                        }
                    }
                }
            }
 
            procMesh.faceZones().clearAddressing();
            procMesh.faceZones().setSize(zoneFaces.size());
            forAll(zoneFaces, zoneI)
            {
                procMesh.faceZones().set
                (
                    zoneI,
                    fz[zoneI].clone
                    (
                        zoneFaces[zoneI].shrink(),          // addressing
                        zoneFaceFlips[zoneI].shrink(),      // flipmap
                        zoneI,
                        procMesh.faceZones()
                    )
                );
            }
 
            if (fz.size())
            {
                // Force writing on all processors
                procMesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
            }
        }
 
        // Cell zones
        {
            const cellZoneMesh& cz = cellZones();
 
            // Go through all the zoned cells and find out if they
            // belong to a zone.  If so, add it to the zone as
            // necessary
            List<DynamicList<label> > zoneCells(cz.size());
 
            // Estimate size
            forAll(zoneCells, zoneI)
            {
                zoneCells[zoneI].setCapacity(cz[zoneI].size() / nProcs_);
            }
 
            forAll(curCellLabels, celli)
            {
                label curCellI = curCellLabels[celli];
 
                label zoneI = cellToZone[curCellI];
 
                if (zoneI >= 0)
                {
                    // Single zone.
                    zoneCells[zoneI].append(celli);
                }
                else if (zoneI == -2)
                {
                    // Multiple zones. Lookup.
                    forAll(cz, zoneI)
                    {
                        label index = cz[zoneI].whichCell(curCellI);
 
                        if (index != -1)
                        {
                            zoneCells[zoneI].append(celli);
                        }
                    }
                }
            }
 
            procMesh.cellZones().clearAddressing();
            procMesh.cellZones().setSize(zoneCells.size());
            forAll(zoneCells, zoneI)
            {
                procMesh.cellZones().set
                (
                    zoneI,
                    cz[zoneI].clone
                    (
                        zoneCells[zoneI].shrink(),
                        zoneI,
                        procMesh.cellZones()
                    )
                );
            }
 
            if (cz.size())
            {
                // Force writing on all processors
                procMesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
            }
        }
 
        // Set the precision of the points data to be min 10
        IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
 
        procMesh.write();
 
        // Write points if pointsInstance differing from facesInstance
        if (facesInstancePointsPtr_.valid())
        {
            pointIOField pointsInstancePoints
            (
                IOobject
                (
                    "points",
                    pointsInstance(),
                    polyMesh::meshSubDir,
                    procMesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE,
                    false
                ),
                xferMove(procPoints)
            );
            pointsInstancePoints.write();
        }
 
 
        // Decompose any sets
        if (decomposeSets)
        {
            forAll(cellSets, i)
            {
                const cellSet& cs = cellSets[i];
                cellSet set(procMesh, cs.name(), cs.size()/nProcs_);
                forAll(curCellLabels, i)
                {
                    if (cs.found(curCellLabels[i]))
                    {
                        set.insert(i);
                    }
                }
                set.write();
            }
            forAll(faceSets, i)
            {
                const faceSet& cs = faceSets[i];
                faceSet set(procMesh, cs.name(), cs.size()/nProcs_);
                forAll(curFaceLabels, i)
                {
                    if (cs.found(mag(curFaceLabels[i])-1))
                    {
                        set.insert(i);
                    }
                }
                set.write();
            }
            forAll(pointSets, i)
            {
                const pointSet& cs = pointSets[i];
                pointSet set(procMesh, cs.name(), cs.size()/nProcs_);
                forAll(curPointLabels, i)
                {
                    if (cs.found(curPointLabels[i]))
                    {
                        set.insert(i);
                    }
                }
                set.write();
            }
        }
 
 
        // Optional hexRef8 data
        hexRef8Data
        (
            IOobject
            (
                "dummy",
                facesInstance(),
                polyMesh::meshSubDir,
                procMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            baseMeshData,
            procCellAddressing_[procI],
            procPointAddressing_[procI]
        ).write();
 
 
        // Statistics
 
        Info<< endl
            << "Processor " << procI << nl
            << "    Number of cells = " << procMesh.nCells()
            << endl;
 
        maxProcCells = max(maxProcCells, procMesh.nCells());
 
        label nBoundaryFaces = 0;
        label nProcPatches = 0;
        label nProcFaces = 0;
 
        forAll(procMesh.boundaryMesh(), patchi)
        {
            if (isA<processorPolyPatch>(procMesh.boundaryMesh()[patchi]))
            {
                const processorPolyPatch& ppp =
                refCast<const processorPolyPatch>
                (
                    procMesh.boundaryMesh()[patchi]
                );
 
                Info<< "    Number of faces shared with processor "
                    << ppp.neighbProcNo() << " = " << ppp.size() << endl;
 
                nProcPatches++;
                nProcFaces += ppp.size();
            }
            else
            {
                nBoundaryFaces += procMesh.boundaryMesh()[patchi].size();
            }
        }
 
        Info<< "    Number of processor patches = " << nProcPatches << nl
            << "    Number of processor faces = " << nProcFaces << nl
            << "    Number of boundary faces = " << nBoundaryFaces << endl;
 
        totProcFaces += nProcFaces;
        totProcPatches += nProcPatches;
        maxProcPatches = max(maxProcPatches, nProcPatches);
        maxProcFaces = max(maxProcFaces, nProcFaces);
 
        // create and write the addressing information
        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();
 
        // Write patch map for backwards compatibility.
        // (= identity map for original patches, -1 for processor patches)
        label nMeshPatches = curPatchSizes.size();
        labelList procBoundaryAddressing(identity(nMeshPatches));
        procBoundaryAddressing.setSize(nMeshPatches+nProcPatches, -1);
 
        labelIOList boundaryProcAddressing
        (
            IOobject
            (
                "boundaryProcAddressing",
                procMesh.facesInstance(),
                procMesh.meshSubDir,
                procMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            procBoundaryAddressing
        );
        boundaryProcAddressing.write();
    }
 
    scalar avgProcCells = scalar(nCells())/nProcs_;
    scalar avgProcPatches = scalar(totProcPatches)/nProcs_;
    scalar avgProcFaces = scalar(totProcFaces)/nProcs_;
 
    // In case of all faces on one processor. Just to avoid division by 0.
    if (totProcPatches == 0)
    {
        avgProcPatches = 1;
    }
    if (totProcFaces == 0)
    {
        avgProcFaces = 1;
    }
 
    Info<< nl
        << "Number of processor faces = " << totProcFaces/2 << nl
        << "Max number of cells = " << maxProcCells
        << " (" << 100.0*(maxProcCells-avgProcCells)/avgProcCells
        << "% above average " << avgProcCells << ")" << nl
        << "Max number of processor patches = " << maxProcPatches
        << " (" << 100.0*(maxProcPatches-avgProcPatches)/avgProcPatches
        << "% above average " << avgProcPatches << ")" << nl
        << "Max number of faces between processors = " << maxProcFaces
        << " (" << 100.0*(maxProcFaces-avgProcFaces)/avgProcFaces
        << "% above average " << avgProcFaces << ")" << nl
        << endl;
 
    return true;
}
 
 
// ************************************************************************* //