inverseDistanceCellCellStencil.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2017-2019 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify i
    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 "inverseDistanceCellCellStencil.H"
#include "addToRunTimeSelectionTable.H"
#include "OBJstream.H"
#include "Time.H"
#include "fvMeshSubset.H"
 
#include "globalIndex.H"
#include "oversetFvPatch.H"
#include "zeroGradientFvPatchFields.H"
#include "syncTools.H"
#include "treeBoundBoxList.H"
#include "waveMethod.H"
 
#include "regionSplit.H"
#include "dynamicOversetFvMesh.H"
//#include "minData.H"
//#include "FaceCellWave.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
namespace cellCellStencils
{
    defineTypeNameAndDebug(inverseDistance, 0);
    addToRunTimeSelectionTable(cellCellStencil, inverseDistance, mesh);
}
}
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
Foam::label Foam::cellCellStencils::inverseDistance::index
(
    const labelVector& nDivs,
    const labelVector& ijk
)
{
    return (ijk[0]*nDivs[1] + ijk[1])*nDivs[2] + ijk[2];
}
 
 
Foam::labelVector Foam::cellCellStencils::inverseDistance::index3
(
    const labelVector& nDivs,
    const label boxI
)
{
    label ij = boxI/nDivs[2];
    label k = boxI-ij*nDivs[2];
    label i = ij/nDivs[1];
    label j = ij-i*nDivs[1];
 
    return labelVector(i, j, k);
}
 
 
Foam::labelVector Foam::cellCellStencils::inverseDistance::index3
(
    const boundBox& bb,
    const labelVector& nDivs,
    const point& pt
)
{
    const vector d(bb.span());
    const point relPt(pt-bb.min());
 
    return labelVector
    (
        floor(relPt[0]/d[0]*nDivs[0]),
        floor(relPt[1]/d[1]*nDivs[1]),
        floor(relPt[2]/d[2]*nDivs[2])
    );
}
 
 
Foam::point Foam::cellCellStencils::inverseDistance::position
(
    const boundBox& bb,
    const labelVector& nDivs,
    const label boxI
)
{
    // Return midpoint of box indicated by boxI
    labelVector ids(index3(nDivs, boxI));
 
    const vector d(bb.span());
    const vector sz(d[0]/nDivs[0], d[1]/nDivs[1], d[2]/nDivs[2]);
 
    return bb.min()+0.5*sz+vector(sz[0]*ids[0], sz[1]*ids[1], sz[2]*ids[2]);
}
 
 
void Foam::cellCellStencils::inverseDistance::fill
(
    PackedList<2>& elems,
    const boundBox& bb,
    const labelVector& nDivs,
    const boundBox& subBb,
    const unsigned int val
)
{
    labelVector minIds(index3(bb, nDivs, subBb.min()));
    labelVector maxIds(index3(bb, nDivs, subBb.max()));
 
    for (direction cmpt = 0; cmpt < 3; cmpt++)
    {
        if (maxIds[cmpt] < 0 || minIds[cmpt] > nDivs[cmpt])
        {
            return;
        }
    }
 
    labelVector maxIndex(labelVector(nDivs[0]-1, nDivs[1]-1, nDivs[2]-1));
    minIds = max(labelVector::zero, minIds);
    maxIds = min(maxIndex, maxIds);
 
    for (label i = minIds[0]; i <= maxIds[0]; i++)
    {
        for (label j = minIds[1]; j <= maxIds[1]; j++)
        {
            for (label k = minIds[2]; k <= maxIds[2]; k++)
            {
                label i1 = index(nDivs, labelVector(i, j, k));
                elems[i1] = val;
            }
        }
    }
}
 
 
void Foam::cellCellStencils::inverseDistance::markBoundaries
(
    const fvMesh& mesh,
    const vector& smallVec,
 
    const boundBox& bb,
    const labelVector& nDivs,
    PackedList<2>& patchTypes,
 
    const labelList& cellMap,
    labelList& patchCellTypes
)
{
    // Mark all voxels that overlap the bounding box of any patch
 
    const fvBoundaryMesh& pbm = mesh.boundary();
 
    patchTypes = patchCellType::OTHER;
 
    // Mark wall boundaries
    forAll(pbm, patchI)
    {
        const fvPatch& fvp = pbm[patchI];
        const labelList& fc = fvp.faceCells();
 
        if (!fvPatch::constraintType(fvp.type()))
        {
            //Info<< "Marking cells on proper patch " << fvp.name()
            //    << " with type " << fvp.type() << endl;
            const polyPatch& pp = fvp.patch();
            forAll(pp, i)
            {
                // Mark in overall patch types
                patchCellTypes[cellMap[fc[i]]] = patchCellType::PATCH;
 
                // Mark in voxel mesh
                boundBox faceBb(pp.points(), pp[i], false);
                faceBb.min() -= smallVec;
                faceBb.max() += smallVec;
 
                if (bb.overlaps(faceBb))
                {
                    fill(patchTypes, bb, nDivs, faceBb, patchCellType::PATCH);
                }
            }
        }
    }
 
    // Override with overset boundaries
    forAll(pbm, patchI)
    {
        const fvPatch& fvp = pbm[patchI];
        const labelList& fc = fvp.faceCells();
 
        if (isA<oversetFvPatch>(fvp))
        {
            //Info<< "Marking cells on overset patch " << fvp.name() << endl;
            const polyPatch& pp = fvp.patch();
            forAll(pp, i)
            {
                // Mark in overall patch types
                patchCellTypes[cellMap[fc[i]]] = patchCellType::OVERSET;
 
                // Mark in voxel mesh
                boundBox faceBb(pp.points(), pp[i], false);
                faceBb.min() -= smallVec;
                faceBb.max() += smallVec;
 
                if (bb.overlaps(faceBb))
                {
                    fill(patchTypes, bb, nDivs, faceBb, patchCellType::OVERSET);
                }
            }
        }
    }
}
 
 
Foam::treeBoundBox Foam::cellCellStencils::inverseDistance::cellBb
(
    const primitiveMesh& mesh,
    const label celli
)
{
    const cellList& cells = mesh.cells();
    const faceList& faces = mesh.faces();
    const pointField& points = mesh.points();
 
    treeBoundBox bb
    (
        vector(GREAT, GREAT, GREAT),
        vector(-GREAT, -GREAT, -GREAT)
    );
 
    const cell& cFaces = cells[celli];
 
    forAll(cFaces, cFacei)
    {
        const face& f = faces[cFaces[cFacei]];
 
        forAll(f, fp)
        {
            const point& p = points[f[fp]];
 
            bb.min() = min(bb.min(), p);
            bb.max() = max(bb.max(), p);
        }
    }
    return bb;
}
 
 
bool Foam::cellCellStencils::inverseDistance::overlaps
(
    const boundBox& bb,
    const labelVector& nDivs,
    const PackedList<2>& vals,
    const treeBoundBox& subBb,
    const unsigned int val
)
{
    // Checks if subBb overlaps any voxel set to val
 
    labelVector minIds(index3(bb, nDivs, subBb.min()));
    labelVector maxIds(index3(bb, nDivs, subBb.max()));
 
    for (direction cmpt = 0; cmpt < 3; cmpt++)
    {
        if (maxIds[cmpt] < 0 || minIds[cmpt] > nDivs[cmpt])
        {
            return false;
        }
    }
 
    labelVector maxIndex(labelVector(nDivs[0]-1, nDivs[1]-1, nDivs[2]-1));
    minIds = max(labelVector::zero, minIds);
    maxIds = min(maxIndex, maxIds);
 
    for (label i = minIds[0]; i <= maxIds[0]; i++)
    {
        for (label j = minIds[1]; j <= maxIds[1]; j++)
        {
            for (label k = minIds[2]; k <= maxIds[2]; k++)
            {
                label i1 = index(nDivs, labelVector(i, j, k));
                if (vals[i1] == patchCellType::PATCH)
                {
                    return true;
                }
            }
        }
    }
    return false;
}
 
 
void Foam::cellCellStencils::inverseDistance::markPatchesAsHoles
(
    PstreamBuffers& pBufs,
 
    const PtrList<fvMeshSubset>& meshParts,
 
    const List<treeBoundBoxList>& patchBb,
    const List<labelVector>& patchDivisions,
    const PtrList<PackedList<2>>& patchParts,
 
    const label srcI,
    const label tgtI,
    labelList& allCellTypes
) const
{
    const treeBoundBoxList& srcPatchBbs = patchBb[srcI];
    const treeBoundBoxList& tgtPatchBbs = patchBb[tgtI];
    const labelList& tgtCellMap = meshParts[tgtI].cellMap();
 
    // 1. do processor-local src-tgt patch overlap
    {
        const treeBoundBox& srcPatchBb = srcPatchBbs[Pstream::myProcNo()];
        const treeBoundBox& tgtPatchBb = tgtPatchBbs[Pstream::myProcNo()];
 
        if (srcPatchBb.overlaps(tgtPatchBb))
        {
            const PackedList<2>& srcPatchTypes = patchParts[srcI];
            const labelVector& zoneDivs = patchDivisions[srcI];
 
            forAll(tgtCellMap, tgtCelli)
            {
                label celli = tgtCellMap[tgtCelli];
                treeBoundBox cBb(cellBb(mesh_, celli));
                cBb.min() -= smallVec_;
                cBb.max() += smallVec_;
 
                if
                (
                    overlaps
                    (
                        srcPatchBb,
                        zoneDivs,
                        srcPatchTypes,
                        cBb,
                        patchCellType::PATCH
                    )
                )
                {
                    allCellTypes[celli] = HOLE;
                }
            }
        }
    }
 
 
    // 2. Send over srcMesh bits that overlap tgt and do calculation
    pBufs.clear();
    for (label procI = 0; procI < Pstream::nProcs(); procI++)
    {
        if (procI != Pstream::myProcNo())
        {
            const treeBoundBox& srcPatchBb = srcPatchBbs[Pstream::myProcNo()];
            const treeBoundBox& tgtPatchBb = tgtPatchBbs[procI];
 
            if (srcPatchBb.overlaps(tgtPatchBb))
            {
                // Send over complete patch voxel map. Tbd: could
                // subset
                UOPstream os(procI, pBufs);
                os << srcPatchBb << patchDivisions[srcI] << patchParts[srcI];
            }
        }
    }
    pBufs.finishedSends();
    for (label procI = 0; procI < Pstream::nProcs(); procI++)
    {
        if (procI != Pstream::myProcNo())
        {
            //const treeBoundBox& srcBb = srcBbs[procI];
            const treeBoundBox& srcPatchBb = srcPatchBbs[procI];
            const treeBoundBox& tgtPatchBb = tgtPatchBbs[Pstream::myProcNo()];
 
            if (srcPatchBb.overlaps(tgtPatchBb))
            {
                UIPstream is(procI, pBufs);
                {
                    treeBoundBox receivedBb(is);
                    if (srcPatchBb != receivedBb)
                    {
                        FatalErrorInFunction
                            << "proc:" << procI
                            << " srcPatchBb:" << srcPatchBb
                            << " receivedBb:" << receivedBb
                            << exit(FatalError);
                    }
                }
                const labelVector zoneDivs(is);
                const PackedList<2> srcPatchTypes(is);
 
                forAll(tgtCellMap, tgtCelli)
                {
                    label celli = tgtCellMap[tgtCelli];
                    treeBoundBox cBb(cellBb(mesh_, celli));
                    cBb.min() -= smallVec_;
                    cBb.max() += smallVec_;
                    if
                    (
                        overlaps
                        (
                            srcPatchBb,
                            zoneDivs,
                            srcPatchTypes,
                            cBb,
                            patchCellType::PATCH
                        )
                    )
                    {
                        allCellTypes[celli] = HOLE;
                    }
                }
            }
        }
    }
}
 
 
bool Foam::cellCellStencils::inverseDistance::betterDonor
(
    const label destMesh,
    const label currentDonorMesh,
    const label newDonorMesh
) const
{
    // This determines for multiple overlapping meshes which one provides
    // the best donors. Is very basic and only looks at indices of meshes:
    // - 'nearest' mesh index wins, i.e. on mesh 0 it preferentially uses donors
    //   from mesh 1 over mesh 2 (if applicable)
    // - if same 'distance' the highest mesh wins. So on mesh 1 it
    //   preferentially uses donors from mesh 2 over mesh 0. This particular
    //   rule helps to avoid some interpolation loops where mesh 1 uses donors
    //   from mesh 0 (usually the background) but mesh 0 then uses
    //   donors from 1.
 
    if (currentDonorMesh == -1)
    {
        return true;
    }
    else
    {
        const label currentDist = mag(currentDonorMesh-destMesh);
        const label newDist = mag(newDonorMesh-destMesh);
 
        if (newDist < currentDist)
        {
            return true;
        }
        else if (newDist == currentDist && newDonorMesh > currentDonorMesh)
        {
            return true;
        }
        else
        {
            return false;
        }
    }
}
 
 
void Foam::cellCellStencils::inverseDistance::markDonors
(
    const globalIndex& globalCells,
    PstreamBuffers& pBufs,
    const PtrList<fvMeshSubset>& meshParts,
    const List<treeBoundBoxList>& meshBb,
 
    const labelList& allCellTypes,
 
    const label srcI,
    const label tgtI,
    labelListList& allStencil,
    labelList& allDonor
) const
{
    const treeBoundBoxList& srcBbs = meshBb[srcI];
    const treeBoundBoxList& tgtBbs = meshBb[tgtI];
 
    const fvMesh& srcMesh = meshParts[srcI].subMesh();
    const labelList& srcCellMap = meshParts[srcI].cellMap();
    const fvMesh& tgtMesh = meshParts[tgtI].subMesh();
    const pointField& tgtCc = tgtMesh.cellCentres();
    const labelList& tgtCellMap = meshParts[tgtI].cellMap();
 
    // 1. do processor-local src/tgt overlap
    {
        labelList tgtToSrcAddr;
        waveMethod::calculate(tgtMesh, srcMesh, tgtToSrcAddr);
        forAll(tgtCellMap, tgtCelli)
        {
            label srcCelli = tgtToSrcAddr[tgtCelli];
            if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
            {
                label celli = tgtCellMap[tgtCelli];
 
                // TBD: check for multiple donors. Maybe better one? For
                //      now check 'nearer' mesh
                if (betterDonor(tgtI, allDonor[celli], srcI))
                {
                    label globalDonor =
                        globalCells.toGlobal(srcCellMap[srcCelli]);
                    allStencil[celli].setSize(1);
                    allStencil[celli][0] = globalDonor;
                    allDonor[celli] = srcI;
                }
            }
        }
    }
 
 
    // 2. Send over tgtMesh bits that overlap src and do calculation on
    //    srcMesh.
 
 
    // (remote) processors where the tgt overlaps my src
    DynamicList<label> tgtOverlapProcs(Pstream::nProcs());
    // (remote) processors where the src overlaps my tgt
    DynamicList<label> srcOverlapProcs(Pstream::nProcs());
    for (label procI = 0; procI < Pstream::nProcs(); procI++)
    {
        if (procI != Pstream::myProcNo())
        {
            if (tgtBbs[procI].overlaps(srcBbs[Pstream::myProcNo()]))
            {
                tgtOverlapProcs.append(procI);
            }
            if (srcBbs[procI].overlaps(tgtBbs[Pstream::myProcNo()]))
            {
                srcOverlapProcs.append(procI);
            }
        }
    }
 
 
 
    // Indices of tgtcells to send over to each processor
    List<DynamicList<label>> tgtSendCells(Pstream::nProcs());
    forAll(srcOverlapProcs, i)
    {
        label procI = srcOverlapProcs[i];
        tgtSendCells[procI].reserve(tgtMesh.nCells()/srcOverlapProcs.size());
    }
 
 
    forAll(tgtCellMap, tgtCelli)
    {
        label celli = tgtCellMap[tgtCelli];
        if (srcOverlapProcs.size())
        {
            treeBoundBox subBb(cellBb(mesh_, celli));
            subBb.min() -= smallVec_;
            subBb.max() += smallVec_;
 
            forAll(srcOverlapProcs, i)
            {
                label procI = srcOverlapProcs[i];
                if (subBb.overlaps(srcBbs[procI]))
                {
                    tgtSendCells[procI].append(tgtCelli);
                }
            }
        }
    }
 
    // Send target cell centres to overlapping processors
    pBufs.clear();
 
    forAll(srcOverlapProcs, i)
    {
        label procI = srcOverlapProcs[i];
        const labelList& cellIDs = tgtSendCells[procI];
 
        UOPstream os(procI, pBufs);
        os << UIndirectList<point>(tgtCc, cellIDs);
    }
    pBufs.finishedSends();
 
    // Receive bits of target processors; find; send back
    (void)srcMesh.tetBasePtIs();
    forAll(tgtOverlapProcs, i)
    {
        label procI = tgtOverlapProcs[i];
 
        UIPstream is(procI, pBufs);
        pointList samples(is);
 
        labelList donors(samples.size(), -1);
        forAll(samples, sampleI)
        {
            const point& sample = samples[sampleI];
            label srcCelli = srcMesh.findCell(sample, polyMesh::CELL_TETS);
            if (srcCelli != -1 && allCellTypes[srcCellMap[srcCelli]] != HOLE)
            {
                donors[sampleI] = globalCells.toGlobal(srcCellMap[srcCelli]);
            }
        }
 
        // Use same pStreamBuffers to send back.
        UOPstream os(procI, pBufs);
        os << donors;
    }
    pBufs.finishedSends();
 
    forAll(srcOverlapProcs, i)
    {
        label procI = srcOverlapProcs[i];
        const labelList& cellIDs = tgtSendCells[procI];
 
        UIPstream is(procI, pBufs);
        labelList donors(is);
 
        if (donors.size() != cellIDs.size())
        {
            FatalErrorInFunction<< "problem : cellIDs:" << cellIDs.size()
                << " donors:" << donors.size() << abort(FatalError);
        }
 
        forAll(donors, donorI)
        {
            label globalDonor = donors[donorI];
 
            if (globalDonor != -1)
            {
                label celli = tgtCellMap[cellIDs[donorI]];
 
                // TBD: check for multiple donors. Maybe better one? For
                //      now check 'nearer' mesh
                if (betterDonor(tgtI, allDonor[celli], srcI))
                {
                    allStencil[celli].setSize(1);
                    allStencil[celli][0] = globalDonor;
                    allDonor[celli] = srcI;
                }
            }
        }
    }
}
 
 
//void Foam::cellCellStencils::inverseDistance::uncompactedRegionSplit
//(
//    const fvMesh& mesh,
//    const globalIndex& globalFaces,
//    const label nZones,
//    const labelList& zoneID,
//    const labelList& cellTypes,
//    const boolList& isBlockedFace,
//    labelList& cellRegion
//) const
//{
//    // Pass 1: locally seed 2 cells per zone (one unblocked, one blocked).
//    // This avoids excessive numbers of front
//
//    // Field on cells and faces.
//    List<minData> cellData(mesh.nCells());
//    List<minData> faceData(mesh.nFaces());
//
//    // Take over blockedFaces by seeding a negative number
//    // (so is always less than the decomposition)
//
//    forAll(isBlockedFace, facei)
//    {
//        if (isBlockedFace[facei])
//        {
//            faceData[facei] = minData(-2);
//        }
//    }
//
//
//    labelList seedFace(nZones, -1);
//
//    const labelList& owner = mesh.faceOwner();
//    const labelList& neighbour = mesh.faceNeighbour();
//
//    forAll(owner, facei)
//    {
//        label own = owner[facei];
//        if (seedFace[zoneID[own]] == -1)
//        {
//            if (cellTypes[own] != HOLE)
//            {
//                const cell& cFaces = mesh.cells()[own];
//                forAll(cFaces, i)
//                {
//                    if (!isBlockedFace[cFaces[i]])
//                    {
//                        seedFace[zoneID[own]] = cFaces[i];
//                    }
//                }
//            }
//        }
//    }
//    forAll(neighbour, facei)
//    {
//        label nei = neighbour[facei];
//        if (seedFace[zoneID[nei]] == -1)
//        {
//            if (cellTypes[nei] != HOLE)
//            {
//                const cell& cFaces = mesh.cells()[nei];
//                forAll(cFaces, i)
//                {
//                    if (!isBlockedFace[cFaces[i]])
//                    {
//                        seedFace[zoneID[nei]] = cFaces[i];
//                    }
//                }
//            }
//        }
//    }
//
//    DynamicList<label> seedFaces(nZones);
//    DynamicList<minData> seedData(seedFaces.size());
//    forAll(seedFace, zonei)
//    {
//        if (seedFace[zonei] != -1)
//        {
//            seedFaces.append(seedFace[zonei]);
//            seedData.append(minData(globalFaces.toGlobal(seedFace[zonei])));
//        }
//    }
//
//    // Propagate information inwards
//    FaceCellWave<minData> deltaCalc
//    (
//        mesh,
//        List<labelPair>(),
//        false,  // disable walking through cyclicAMI for backwards
//                // compatibility
//        seedFaces,
//        seedData,
//        faceData,
//        cellData,
//        mesh.globalData().nTotalCells()+1
//    );
//
//    // Extract
//    cellRegion.setSize(mesh.nCells());
//    forAll(cellRegion, celli)
//    {
//        if (cellData[celli].valid(deltaCalc.data()))
//        {
//            cellRegion[celli] = cellData[celli].data();
//        }
//        else
//        {
//            // Unvisited cell -> only possible if surrounded by blocked faces.
//            // If so make up region from any of the faces
//            const cell& cFaces = mesh.cells()[celli];
//            label facei = cFaces[0];
//            cellRegion[celli] = globalFaces.toGlobal(facei);
//        }
//    }
//}
//Foam::autoPtr<Foam::globalIndex>
//Foam::cellCellStencils::inverseDistance::compactedRegionSplit
//(
//    const fvMesh& mesh,
//    const globalIndex& globalRegions,
//    labelList& cellRegion
//) const
//{
//    // Now our cellRegion will have
//    // - non-local regions (i.e. originating from other processors)
//    // - non-compact locally originating regions
//    // so we'll need to compact
//
//    // 4a: count per originating processor the number of regions
//    labelList nOriginating(Pstream::nProcs(), Zero);
//    {
//        labelHashSet haveRegion(mesh.nCells()/8);
//
//        forAll(cellRegion, celli)
//        {
//            label region = cellRegion[celli];
//
//            // Count originating processor. Use isLocal as efficiency since
//            // most cells are locally originating.
//            if (globalRegions.isLocal(region))
//            {
//                if (haveRegion.insert(region))
//                {
//                    nOriginating[Pstream::myProcNo()]++;
//                }
//            }
//            else
//            {
//                label proci = globalRegions.whichProcID(region);
//                if (haveRegion.insert(region))
//                {
//                    nOriginating[proci]++;
//                }
//            }
//        }
//    }
//
//    if (debug)
//    {
//        Pout<< "Counted " << nOriginating[Pstream::myProcNo()]
//            << " local regions." << endl;
//    }
//
//
//    // Global numbering for compacted local regions
//    autoPtr<globalIndex> globalCompactPtr
//    (
//        new globalIndex(nOriginating[Pstream::myProcNo()])
//    );
//    const globalIndex& globalCompact = globalCompactPtr();
//
//
//    // 4b: renumber
//    // Renumber into compact indices. Note that since we've already made
//    // all regions global we now need a Map to store the compacting
//    // information
//    // instead of a labelList - otherwise we could have used a straight
//    // labelList.
//
//    // Local compaction map
//    Map<label> globalToCompact(2*nOriginating[Pstream::myProcNo()]);
//    // Remote regions we want the compact number for
//    List<labelHashSet> nonLocal(Pstream::nProcs());
//    forAll(nonLocal, proci)
//    {
//        if (proci != Pstream::myProcNo())
//        {
//            nonLocal[proci].resize(2*nOriginating[proci]);
//        }
//    }
//
//    forAll(cellRegion, celli)
//    {
//        label region = cellRegion[celli];
//        if (globalRegions.isLocal(region))
//        {
//            // Insert new compact region (if not yet present)
//            globalToCompact.insert
//            (
//                region,
//                globalCompact.toGlobal(globalToCompact.size())
//            );
//        }
//        else
//        {
//            nonLocal[globalRegions.whichProcID(region)].insert(region);
//        }
//    }
//
//
//    // Now we have all the local regions compacted. Now we need to get the
//    // non-local ones from the processors to whom they are local.
//    // Convert the nonLocal (labelHashSets) to labelLists.
//
//    labelListList sendNonLocal(Pstream::nProcs());
//    forAll(sendNonLocal, proci)
//    {
//        sendNonLocal[proci] = nonLocal[proci].toc();
//    }
//
//    if (debug)
//    {
//        forAll(sendNonLocal, proci)
//        {
//            Pout<< "    from processor " << proci
//                << " want " << sendNonLocal[proci].size()
//                << " region numbers." << endl;
//        }
//        Pout<< endl;
//    }
//
//
//    // Get the wanted region labels into recvNonLocal
//    labelListList recvNonLocal;
//    Pstream::exchange<labelList, label>(sendNonLocal, recvNonLocal);
//
//    // Now we have the wanted compact region labels that proci wants in
//    // recvNonLocal[proci]. Construct corresponding list of compact
//    // region labels to send back.
//
//    labelListList sendWantedLocal(Pstream::nProcs());
//    forAll(recvNonLocal, proci)
//    {
//        const labelList& nonLocal = recvNonLocal[proci];
//        sendWantedLocal[proci].setSize(nonLocal.size());
//
//        forAll(nonLocal, i)
//        {
//            sendWantedLocal[proci][i] = globalToCompact[nonLocal[i]];
//        }
//    }
//
//
//    // Send back (into recvNonLocal)
//    recvNonLocal.clear();
//    Pstream::exchange<labelList, label>(sendWantedLocal, recvNonLocal);
//    sendWantedLocal.clear();
//
//    // Now recvNonLocal contains for every element in setNonLocal the
//    // corresponding compact number. Insert these into the local compaction
//    // map.
//
//    forAll(recvNonLocal, proci)
//    {
//        const labelList& wantedRegions = sendNonLocal[proci];
//        const labelList& compactRegions = recvNonLocal[proci];
//
//        forAll(wantedRegions, i)
//        {
//            globalToCompact.insert(wantedRegions[i], compactRegions[i]);
//        }
//    }
//
//    // Finally renumber the regions
//    forAll(cellRegion, celli)
//    {
//        cellRegion[celli] = globalToCompact[cellRegion[celli]];
//    }
//
//    return globalCompactPtr;
//}
 
 
void Foam::cellCellStencils::inverseDistance::findHoles
(
    const globalIndex& globalCells,
    const fvMesh& mesh,
    const labelList& zoneID,
    const labelListList& stencil,
    labelList& cellTypes
) const
{
    const fvBoundaryMesh& pbm = mesh.boundary();
    const labelList& own = mesh.faceOwner();
    const labelList& nei = mesh.faceNeighbour();
 
 
    // The input cellTypes will be
    // - HOLE           : cell part covered by other-mesh patch
    // - INTERPOLATED   : cell fully covered by other-mesh patch
    //                    or next to 'overset' patch
    // - CALCULATED     : otherwise
    //
    // so we start a walk from our patches and any cell we cannot reach
    // (because we walk is stopped by other-mesh patch) is a hole.
 
 
    DebugInfo<< FUNCTION_NAME << " : Starting hole flood filling" << endl;
 
    DebugInfo<< FUNCTION_NAME << " : Starting hole cells : "
        << findIndices(cellTypes, HOLE).size() << endl;
 
    boolList isBlockedFace(mesh.nFaces(), false);
    label nBlocked = 0;
 
    for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
    {
        label ownType = cellTypes[own[faceI]];
        label neiType = cellTypes[nei[faceI]];
        if
        (
             (ownType == HOLE && neiType != HOLE)
          || (ownType != HOLE && neiType == HOLE)
        )
        {
            isBlockedFace[faceI] = true;
            nBlocked++;
        }
    }
    DebugInfo<< FUNCTION_NAME << " : Marked internal hole boundaries : "
        << nBlocked << endl;
 
 
    labelList nbrCellTypes;
    syncTools::swapBoundaryCellList(mesh, cellTypes, nbrCellTypes);
 
    for (label faceI = mesh.nInternalFaces(); faceI < mesh.nFaces(); faceI++)
    {
        label ownType = cellTypes[own[faceI]];
        label neiType = nbrCellTypes[faceI-mesh.nInternalFaces()];
 
        if
        (
             (ownType == HOLE && neiType != HOLE)
          || (ownType != HOLE && neiType == HOLE)
        )
        {
            isBlockedFace[faceI] = true;
            nBlocked++;
        }
    }
 
    DebugInfo<< FUNCTION_NAME << " : Marked all hole boundaries : "
        << nBlocked << endl;
 
    // Determine regions
    regionSplit cellRegion(mesh, isBlockedFace);
    const label nRegions = cellRegion.nRegions();
 
    //labelList cellRegion;
    //label nRegions = -1;
    //{
    //    const globalIndex globalFaces(mesh.nFaces());
    //    uncompactedRegionSplit
    //    (
    //        mesh,
    //        globalFaces,
    //        gMax(zoneID)+1,
    //        zoneID,
    //        cellTypes,
    //        isBlockedFace,
    //        cellRegion
    //    );
    //    autoPtr<globalIndex> globalRegions
    //    (
    //        compactedRegionSplit
    //        (
    //            mesh,
    //            globalFaces,
    //            cellRegion
    //        )
    //    );
    //    nRegions = globalRegions().size();
    //}
    DebugInfo<< FUNCTION_NAME << " : Determined regions : "
        << nRegions << endl;
 
    //Info<< typeName << " : detected " << nRegions
    //    << " mesh regions after overset" << nl << endl;
 
 
 
    // Now we'll have a mesh split according to where there are cells
    // covered by the other-side patches. See what we can reach from our
    // real patches
 
    //  0 : region not yet determined
    //  1 : borders blockage so is not ok (but can be overridden by real
    //      patch)
    //  2 : has real patch in it so is reachable
    labelList regionType(nRegions, Zero);
 
 
    // See if any regions borders blockage. Note: isBlockedFace is already
    // parallel synchronised.
    {
        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
        {
            if (isBlockedFace[faceI])
            {
                label ownRegion = cellRegion[own[faceI]];
 
                if (cellTypes[own[faceI]] != HOLE)
                {
                    if (regionType[ownRegion] == 0)
                    {
                        regionType[ownRegion] = 1;
                    }
                }
 
                label neiRegion = cellRegion[nei[faceI]];
 
                if (cellTypes[nei[faceI]] != HOLE)
                {
                    if (regionType[neiRegion] == 0)
                    {
                        regionType[neiRegion] = 1;
                    }
                }
            }
        }
        for
        (
            label faceI = mesh.nInternalFaces();
            faceI < mesh.nFaces();
            faceI++
        )
        {
            if (isBlockedFace[faceI])
            {
                label ownRegion = cellRegion[own[faceI]];
 
                if (regionType[ownRegion] == 0)
                {
                    regionType[ownRegion] = 1;
                }
            }
        }
    }
 
 
    // Override with real patches
    forAll(pbm, patchI)
    {
        const fvPatch& fvp = pbm[patchI];
 
        if (isA<oversetFvPatch>(fvp))
        {}
        else if (!fvPatch::constraintType(fvp.type()))
        {
            const labelList& fc = fvp.faceCells();
            forAll(fc, i)
            {
                label regionI = cellRegion[fc[i]];
 
                if (cellTypes[fc[i]] != HOLE && regionType[regionI] != 2)
                {
                    regionType[regionI] = 2;
                }
            }
        }
    }
 
    DebugInfo<< FUNCTION_NAME << " : Done local analysis" << endl;
 
    // Now we've handled
    // - cells next to blocked cells
    // - coupled boundaries
    // Only thing to handle is the interpolation between regions
 
 
    labelListList compactStencil(stencil);
    List<Map<label>> compactMap;
    mapDistribute map(globalCells, compactStencil, compactMap);
 
    DebugInfo<< FUNCTION_NAME << " : Converted stencil into compact form"
        << endl;
 
 
    while (true)
    {
        // Synchronise region status on processors
        // (could instead swap status through processor patches)
        Pstream::listCombineGather(regionType, maxEqOp<label>());
        Pstream::listCombineScatter(regionType);
 
        DebugInfo<< FUNCTION_NAME << " : Gathered region type" << endl;
 
        // Communicate region status through interpolative cells
        labelList cellRegionType(labelUIndList(regionType, cellRegion));
        map.distribute(cellRegionType);
 
        DebugInfo<< FUNCTION_NAME << " : Interpolated region type" << endl;
 
 
 
        label nChanged = 0;
        forAll(pbm, patchI)
        {
            const fvPatch& fvp = pbm[patchI];
 
            if (isA<oversetFvPatch>(fvp))
            {
                const labelUList& fc = fvp.faceCells();
                forAll(fc, i)
                {
                    label cellI = fc[i];
                    label regionI = cellRegion[cellI];
 
                    if (regionType[regionI] != 2)
                    {
                        const labelList& slots = compactStencil[cellI];
                        forAll(slots, i)
                        {
                            label otherType = cellRegionType[slots[i]];
 
                            if (otherType == 2)
                            {
                                //Pout<< "Reachable through interpolation : "
                                //    << regionI << " at cell "
                                //    << mesh.cellCentres()[cellI] << endl;
                                regionType[regionI] = 2;
                                nChanged++;
                                break;
                            }
                        }
                    }
                }
            }
        }
 
        reduce(nChanged, sumOp<label>());
        DebugInfo<< FUNCTION_NAME << " : Determined regions changed : "
            << nChanged << endl;
 
        if (nChanged == 0)
        {
            break;
        }
    }
 
 
    // See which regions have not been visited (regionType == 1)
    forAll(cellRegion, cellI)
    {
        label type = regionType[cellRegion[cellI]];
        if (type == 1 && cellTypes[cellI] != HOLE)
        {
            cellTypes[cellI] = HOLE;
        }
    }
    DebugInfo<< FUNCTION_NAME << " : Finished hole flood filling" << endl;
}
 
 
void Foam::cellCellStencils::inverseDistance::seedCell
(
    const label cellI,
    const scalar wantedFraction,
    bitSet& isFront,
    scalarField& fraction
) const
{
    const cell& cFaces = mesh_.cells()[cellI];
    forAll(cFaces, i)
    {
        label nbrFacei = cFaces[i];
        if (fraction[nbrFacei] < wantedFraction)
        {
            fraction[nbrFacei] = wantedFraction;
            isFront.set(nbrFacei);
        }
    }
}
 
 
void Foam::cellCellStencils::inverseDistance::walkFront
(
    const scalar layerRelax,
    const labelListList& allStencil,
    labelList& allCellTypes,
    scalarField& allWeight
) const
{
    // Current front
    bitSet isFront(mesh_.nFaces());
 
    const fvBoundaryMesh& fvm = mesh_.boundary();
 
 
    // 'overset' patches
 
    forAll(fvm, patchI)
    {
        if (isA<oversetFvPatch>(fvm[patchI]))
        {
            const labelList& fc = fvm[patchI].faceCells();
            forAll(fc, i)
            {
                label cellI = fc[i];
                if (allCellTypes[cellI] == INTERPOLATED)
                {
                    // Note that acceptors might have been marked hole if
                    // there are no donors in which case we do not want to
                    // walk this out. This is an extreme situation.
                    isFront.set(fvm[patchI].start()+i);
                }
            }
        }
    }
 
 
    // Outside of 'hole' region
    {
        const labelList& own = mesh_.faceOwner();
        const labelList& nei = mesh_.faceNeighbour();
 
        for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
        {
            label ownType = allCellTypes[own[faceI]];
            label neiType = allCellTypes[nei[faceI]];
            if
            (
                 (ownType == HOLE && neiType != HOLE)
              || (ownType != HOLE && neiType == HOLE)
            )
            {
                //Pout<< "Front at face:" << faceI
                //    << " at:" << mesh_.faceCentres()[faceI] << endl;
                isFront.set(faceI);
            }
        }
 
        labelList nbrCellTypes;
        syncTools::swapBoundaryCellList(mesh_, allCellTypes, nbrCellTypes);
 
        for
        (
            label faceI = mesh_.nInternalFaces();
            faceI < mesh_.nFaces();
            faceI++
        )
        {
            label ownType = allCellTypes[own[faceI]];
            label neiType = nbrCellTypes[faceI-mesh_.nInternalFaces()];
 
            if
            (
                 (ownType == HOLE && neiType != HOLE)
              || (ownType != HOLE && neiType == HOLE)
            )
            {
                //Pout<< "Front at coupled face:" << faceI
                //    << " at:" << mesh_.faceCentres()[faceI] << endl;
                isFront.set(faceI);
            }
        }
    }
 
 
    // Current interpolation fraction
    scalarField fraction(mesh_.nFaces(), Zero);
 
    forAll(isFront, faceI)
    {
        if (isFront.test(faceI))
        {
            fraction[faceI] = 1.0;
        }
    }
 
 
    while (returnReduce(isFront.any(), orOp<bool>()))
    {
        // Interpolate cells on front
        bitSet newIsFront(mesh_.nFaces());
        scalarField newFraction(fraction);
        forAll(isFront, faceI)
        {
            if (isFront.test(faceI))
            {
                label own = mesh_.faceOwner()[faceI];
                if (allCellTypes[own] != HOLE)
                {
                    if (allWeight[own] < fraction[faceI])
                    {
                        // Cell wants to become interpolated (if sufficient
                        // stencil, otherwise becomes hole)
                        if (allStencil[own].size())
                        {
                            allWeight[own] = fraction[faceI];
                            allCellTypes[own] = INTERPOLATED;
                            // Add faces of cell (with lower weight) as new
                            // front
                            seedCell
                            (
                                own,
                                fraction[faceI]-layerRelax,
                                newIsFront,
                                newFraction
                            );
                        }
                        else
                        {
                            allWeight[own] = 0.0;
                            allCellTypes[own] = HOLE;
                            // Add faces of cell as new front
                            seedCell
                            (
                                own,
                                1.0,
                                newIsFront,
                                newFraction
                            );
                        }
                    }
                }
                if (mesh_.isInternalFace(faceI))
                {
                    label nei = mesh_.faceNeighbour()[faceI];
                    if (allCellTypes[nei] != HOLE)
                    {
                        if (allWeight[nei] < fraction[faceI])
                        {
                            if (allStencil[nei].size())
                            {
                                allWeight[nei] = fraction[faceI];
                                allCellTypes[nei] = INTERPOLATED;
                                seedCell
                                (
                                    nei,
                                    fraction[faceI]-layerRelax,
                                    newIsFront,
                                    newFraction
                                );
                            }
                            else
                            {
                                allWeight[nei] = 0.0;
                                allCellTypes[nei] = HOLE;
                                seedCell
                                (
                                    nei,
                                    1.0,
                                    newIsFront,
                                    newFraction
                                );
                            }
                        }
                    }
                }
            }
        }
 
        syncTools::syncFaceList(mesh_, newIsFront, orEqOp<unsigned int>());
        syncTools::syncFaceList(mesh_, newFraction, maxEqOp<scalar>());
 
        isFront.transfer(newIsFront);
        fraction.transfer(newFraction);
    }
}
 
 
void Foam::cellCellStencils::inverseDistance::stencilWeights
(
    const point& sample,
    const pointList& donorCcs,
    scalarList& weights
) const
{
    // Inverse-distance weighting
 
    weights.setSize(donorCcs.size());
    scalar sum = 0.0;
    forAll(donorCcs, i)
    {
        scalar d = mag(sample-donorCcs[i]);
 
        if (d > ROOTVSMALL)
        {
            weights[i] = 1.0/d;
            sum += weights[i];
        }
        else
        {
            // Short circuit
            weights = 0.0;
            weights[i] = 1.0;
            return;
        }
    }
    forAll(weights, i)
    {
        weights[i] /= sum;
    }
}
 
 
void Foam::cellCellStencils::inverseDistance::createStencil
(
    const globalIndex& globalCells
)
{
    // Send cell centre back to donor
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    // The complication is that multiple acceptors need the same donor
    // (but with different weights obviously)
    // So we do multi-pass:
    // - send over cc of acceptor for which we want stencil.
    //   Consistently choose the acceptor with smallest magSqr in case of
    //   multiple acceptors for the containing cell/donor.
    // - find the cell-cells and weights for the donor
    // - send back together with the acceptor cc
    // - use the acceptor cc to see if it was 'me' that sent it. If so
    //   mark me as complete so it doesn't get involved in the next loop.
    // - loop until all resolved.
 
    // Special value for unused points
    const vector greatPoint(GREAT, GREAT, GREAT);
 
    boolList isValidDonor(mesh_.nCells(), true);
    forAll(cellTypes_, celli)
    {
        if (cellTypes_[celli] == HOLE)
        {
            isValidDonor[celli] = false;
        }
    }
 
 
    // Has acceptor been handled already?
    bitSet doneAcceptor(interpolationCells_.size());
 
    while (true)
    {
        pointField samples(cellInterpolationMap().constructSize(), greatPoint);
 
        // Fill remote slots (override old content). We'll find out later
        // on which one has won and mark this one in doneAcceptor.
        label nSamples = 0;
        forAll(interpolationCells_, i)
        {
            if (!doneAcceptor[i])
            {
                label cellI = interpolationCells_[i];
                const point& cc = mesh_.cellCentres()[cellI];
                const labelList& slots = cellStencil_[cellI];
 
                if (slots.size() != 1)
                {
                    FatalErrorInFunction<< "Problem:" << slots
                        << abort(FatalError);
                }
 
                forAll(slots, slotI)
                {
                    label elemI = slots[slotI];
                    //Pout<< "    acceptor:" << cellI
                    //    << " at:" << mesh_.cellCentres()[cellI]
                    //    << " global:" << globalCells.toGlobal(cellI)
                    //    << " foudn in donor:" << elemI << endl;
                    minMagSqrEqOp<point>()(samples[elemI], cc);
                }
                nSamples++;
            }
        }
 
 
        if (returnReduce(nSamples, sumOp<label>()) == 0)
        {
            break;
        }
 
        // Send back to donor. Make sure valid point takes priority
        mapDistributeBase::distribute<point, minMagSqrEqOp<point>, flipOp>
        (
            Pstream::commsTypes::nonBlocking,
            List<labelPair>(),
            mesh_.nCells(),
            cellInterpolationMap().constructMap(),
            false,
            cellInterpolationMap().subMap(),
            false,
            samples,
            minMagSqrEqOp<point>(),
            flipOp(),                               // negateOp
            greatPoint                              // nullValue
        );
 
        // All the donor cells will now have a valid cell centre. Construct a
        // stencil for these.
 
        DynamicList<label> donorCells(mesh_.nCells());
        forAll(samples, cellI)
        {
            if (samples[cellI] != greatPoint)
            {
                donorCells.append(cellI);
            }
        }
 
 
        // Get neighbours (global cell and centre) of donorCells.
        labelListList donorCellCells(mesh_.nCells());
        pointListList donorCellCentres(mesh_.nCells());
        globalCellCells
        (
            globalCells,
            mesh_,
            isValidDonor,
            donorCells,
            donorCellCells,
            donorCellCentres
        );
 
        // Determine the weights.
        scalarListList donorWeights(mesh_.nCells());
        forAll(donorCells, i)
        {
            label cellI = donorCells[i];
            const pointList& donorCentres = donorCellCentres[cellI];
            stencilWeights
            (
                samples[cellI],
                donorCentres,
                donorWeights[cellI]
            );
        }
 
        // Transfer the information back to the acceptor:
        // - donorCellCells : stencil (with first element the original donor)
        // - donorWeights : weights for donorCellCells
        cellInterpolationMap().distribute(donorCellCells);
        cellInterpolationMap().distribute(donorWeights);
        cellInterpolationMap().distribute(samples);
 
        // Check which acceptor has won and transfer
        forAll(interpolationCells_, i)
        {
            if (!doneAcceptor[i])
            {
                label cellI = interpolationCells_[i];
                const labelList& slots = cellStencil_[cellI];
 
                if (slots.size() != 1)
                {
                    FatalErrorInFunction << "Problem:" << slots
                        << abort(FatalError);
                }
 
                label slotI = slots[0];
 
                // Important: check if the stencil is actually for this cell
                if (samples[slotI] == mesh_.cellCentres()[cellI])
                {
                    cellStencil_[cellI].transfer(donorCellCells[slotI]);
                    cellInterpolationWeights_[cellI].transfer
                    (
                        donorWeights[slotI]
                    );
                    // Mark cell as being done so it does not get sent over
                    // again.
                    doneAcceptor.set(i);
                }
            }
        }
    }
 
    // Re-do the mapDistribute
    List<Map<label>> compactMap;
    cellInterpolationMap_.reset
    (
        new mapDistribute
        (
            globalCells,
            cellStencil_,
            compactMap
        )
    );
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::cellCellStencils::inverseDistance::inverseDistance
(
    const fvMesh& mesh,
    const dictionary& dict,
    const bool doUpdate
)
:
    cellCellStencil(mesh),
    dict_(dict),
    smallVec_(Zero),
    cellTypes_(labelList(mesh.nCells(), CALCULATED)),
    interpolationCells_(0),
    cellInterpolationMap_(),
    cellStencil_(0),
    cellInterpolationWeights_(0),
    cellInterpolationWeight_
    (
        IOobject
        (
            "cellInterpolationWeight",
            mesh_.facesInstance(),
            mesh_,
            IOobject::NO_READ,
            IOobject::NO_WRITE,
            false
        ),
        mesh_,
        dimensionedScalar(dimless, Zero),
        zeroGradientFvPatchScalarField::typeName
    )
{
    // Protect local fields from interpolation
    nonInterpolatedFields_.insert("cellInterpolationWeight");
    nonInterpolatedFields_.insert("cellTypes");
    nonInterpolatedFields_.insert("maxMagWeight");
 
    // For convenience also suppress frequently used displacement field
    nonInterpolatedFields_.insert("cellDisplacement");
    nonInterpolatedFields_.insert("grad(cellDisplacement)");
    const word w("snGradCorr(cellDisplacement)");
    const word d("((viscosity*faceDiffusivity)*magSf)");
    nonInterpolatedFields_.insert("surfaceIntegrate(("+d+"*"+w+"))");
 
    // Read zoneID
    this->zoneID();
 
    // Read old-time cellTypes
    IOobject io
    (
        "cellTypes",
        mesh_.time().timeName(),
        mesh_,
        IOobject::READ_IF_PRESENT,
        IOobject::NO_WRITE,
        false
    );
    if (io.typeHeaderOk<volScalarField>(true))
    {
        if (debug)
        {
            Pout<< "Reading cellTypes from time " << mesh_.time().timeName()
                << endl;
        }
 
        const volScalarField volCellTypes(io, mesh_);
        forAll(volCellTypes, celli)
        {
            // Round to integer
            cellTypes_[celli] = volCellTypes[celli];
        }
    }
 
    if (doUpdate)
    {
        update();
    }
}
 
 
// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
Foam::cellCellStencils::inverseDistance::~inverseDistance()
{}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
bool Foam::cellCellStencils::inverseDistance::update()
{
    scalar layerRelax(dict_.lookupOrDefault("layerRelax", 1.0));
 
    scalar tol = dict_.lookupOrDefault("tolerance", 1e-10);
    smallVec_ = mesh_.bounds().span()*tol;
 
    const labelIOList& zoneID = this->zoneID();
 
    label nZones = gMax(zoneID)+1;
    labelList nCellsPerZone(nZones, Zero);
    forAll(zoneID, cellI)
    {
        nCellsPerZone[zoneID[cellI]]++;
    }
    Pstream::listCombineGather(nCellsPerZone, plusEqOp<label>());
    Pstream::listCombineScatter(nCellsPerZone);
 
 
    const boundBox& allBb(mesh_.bounds());
 
 
    PtrList<fvMeshSubset> meshParts(nZones);
    List<treeBoundBoxList> meshBb(nZones);
 
    // Determine zone meshes and bounding boxes
    {
        // Per processor, per zone the bounding box
        List<treeBoundBoxList> procBb(Pstream::nProcs());
        procBb[Pstream::myProcNo()].setSize(nZones);
 
        forAll(meshParts, zonei)
        {
            meshParts.set
            (
                zonei,
                new fvMeshSubset(mesh_, zonei, zoneID)
            );
            const fvMesh& subMesh = meshParts[zonei].subMesh();
 
            // Trigger early evaluation of mesh dimension (in case there are
            // zero cells in mesh)
            (void)subMesh.nGeometricD();
 
            if (subMesh.nPoints())
            {
                procBb[Pstream::myProcNo()][zonei] =
                    treeBoundBox(subMesh.points());
                procBb[Pstream::myProcNo()][zonei].inflate(1e-6);
            }
            else
            {
                // No part of zone on this processor. Make up bb.
                procBb[Pstream::myProcNo()][zonei] = treeBoundBox
                (
                    allBb.min() - 2*allBb.span(),
                    allBb.min() - allBb.span()
                );
                procBb[Pstream::myProcNo()][zonei].inflate(1e-6);
            }
        }
 
        Pstream::gatherList(procBb);
        Pstream::scatterList(procBb);
 
        // Move local bounding boxes to per-mesh indexing
        forAll(meshBb, zoneI)
        {
            treeBoundBoxList& bbs = meshBb[zoneI];
            bbs.setSize(Pstream::nProcs());
            forAll(procBb, procI)
            {
                bbs[procI] = procBb[procI][zoneI];
            }
        }
    }
 
 
    // Determine patch bounding boxes. These are either global and provided
    // by the user or processor-local as a copy of the mesh bounding box.
 
    List<treeBoundBoxList> patchBb(nZones);
    List<labelVector> patchDivisions(nZones);
    PtrList<PackedList<2>> patchParts(nZones);
    labelList allPatchTypes(mesh_.nCells(), OTHER);
 
    {
        treeBoundBox globalPatchBb;
        if (dict_.readIfPresent("searchBox", globalPatchBb))
        {
            // All processors, all zones have the same bounding box
            patchBb = treeBoundBoxList(Pstream::nProcs(), globalPatchBb);
        }
        else
        {
            // Use the meshBb (differing per zone, per processor)
            patchBb = meshBb;
        }
    }
 
    {
        labelVector globalDivs;
        if (dict_.readIfPresent("searchBoxDivisions", globalDivs))
        {
            patchDivisions = globalDivs;
        }
        else
        {
            const labelVector& dim = mesh_.geometricD();
            label nDivs = -1;
            if (mesh_.nGeometricD() == 1)
            {
                nDivs = mesh_.nCells();
            }
            else if (mesh_.nGeometricD() == 2)
            {
                nDivs = label(Foam::sqrt(scalar(mesh_.nCells())));
            }
            else
            {
                nDivs = label(Foam::cbrt(scalar(mesh_.nCells())));
            }
 
            labelVector v(nDivs, nDivs, nDivs);
            forAll(dim, i)
            {
                if (dim[i] == -1)
                {
                    v[i] = 1;
                }
            }
            patchDivisions = v;
        }
    }
 
    forAll(patchParts, zoneI)
    {
        patchParts.set
        (
            zoneI,
            new PackedList<2>
            (
                patchDivisions[zoneI][0]
               *patchDivisions[zoneI][1]
               *patchDivisions[zoneI][2]
            )
        );
        markBoundaries
        (
            meshParts[zoneI].subMesh(),
            smallVec_,
 
            patchBb[zoneI][Pstream::myProcNo()],
            patchDivisions[zoneI],
            patchParts[zoneI],
 
            meshParts[zoneI].cellMap(),
            allPatchTypes
        );
    }
 
 
    // Print a bit
    {
        Info<< type() << " : detected " << nZones
            << " mesh regions" << endl;
        Info<< incrIndent;
        forAll(nCellsPerZone, zoneI)
        {
            Info<< indent<< "zone:" << zoneI
                << " nCells:" << nCellsPerZone[zoneI]
                << "  voxels:" << patchDivisions[zoneI]
                << " bb:" << patchBb[zoneI][Pstream::myProcNo()]
                << endl;
        }
        Info<< decrIndent;
    }
 
 
    // Current best guess for cells. Includes best stencil. Weights should
    // add up to volume.
    labelList allCellTypes(mesh_.nCells(), CALCULATED);
    labelListList allStencil(mesh_.nCells());
    // zoneID of donor
    labelList allDonorID(mesh_.nCells(), -1);
 
    const globalIndex globalCells(mesh_.nCells());
 
    PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
 
    // Mark holes (in allCellTypes)
    for (label srcI = 0; srcI < meshParts.size()-1; srcI++)
    {
        for (label tgtI = srcI+1; tgtI < meshParts.size(); tgtI++)
        {
            markPatchesAsHoles
            (
                pBufs,
 
                meshParts,
 
                patchBb,
                patchDivisions,
                patchParts,
 
                srcI,
                tgtI,
                allCellTypes
            );
            markPatchesAsHoles
            (
                pBufs,
 
                meshParts,
 
                patchBb,
                patchDivisions,
                patchParts,
 
                tgtI,
                srcI,
                allCellTypes
            );
        }
    }
 
    // Find donors (which are not holes) in allStencil, allDonorID
    for (label srcI = 0; srcI < meshParts.size()-1; srcI++)
    {
        for (label tgtI = srcI+1; tgtI < meshParts.size(); tgtI++)
        {
            markDonors
            (
                globalCells,
                pBufs,
                meshParts,
                meshBb,
                allCellTypes,
 
                tgtI,
                srcI,
                allStencil,
                allDonorID
            );
            markDonors
            (
                globalCells,
                pBufs,
                meshParts,
                meshBb,
                allCellTypes,
 
                srcI,
                tgtI,
                allStencil,
                allDonorID
            );
        }
    }
 
    if (debug)
    {
        tmp<volScalarField> tfld
        (
            createField(mesh_, "allCellTypes", allCellTypes)
        );
        tfld().write();
    }
 
    // Use the patch types and weights to decide what to do
    forAll(allPatchTypes, cellI)
    {
        if (allCellTypes[cellI] != HOLE)
        {
            switch (allPatchTypes[cellI])
            {
                case OVERSET:
                {
                    // Require interpolation. See if possible.
 
                    if (allStencil[cellI].size())
                    {
                        allCellTypes[cellI] = INTERPOLATED;
                    }
                    else
                    {
                        // If there are no donors we can either set the
                        // acceptor cell to 'hole' i.e. nothing gets calculated
                        // if the acceptor cells go outside the donor meshes or
                        // we can set it to 'calculated' to have something
                        // like an acmi type behaviour where only covered
                        // acceptor cell use the interpolation and non-covered
                        // become calculated. Uncomment below line. Note: this
                        // should be switchable?
                        //allCellTypes[cellI] = CALCULATED;
 
                        allCellTypes[cellI] = HOLE;
                    }
                }
            }
        }
    }
 
    if (debug)
    {
        tmp<volScalarField> tfld
        (
            createField(mesh_, "allCellTypes_patch", allCellTypes)
        );
        tfld().write();
    }
 
 
    // Check previous iteration cellTypes_ for any hole->calculated changes
    // If so set the cell either to interpolated (if there are donors) or
    // holes (if there are no donors). Note that any interpolated cell might
    // still be overwritten by the flood filling
    {
        label nCalculated = 0;
 
        forAll(cellTypes_, celli)
        {
            if (allCellTypes[celli] == CALCULATED && cellTypes_[celli] == HOLE)
            {
                if (allStencil[celli].size() == 0)
                {
                    // Reset to hole
                    allCellTypes[celli] = HOLE;
                    allStencil[celli].clear();
                }
                else
                {
                    allCellTypes[celli] = INTERPOLATED;
                    nCalculated++;
                }
            }
        }
 
        if (debug)
        {
            Pout<< "Detected " << nCalculated << " cells changing from hole"
                << " to calculated. Changed to interpolated"
                << endl;
        }
    }
 
 
    // Mark unreachable bits
    findHoles(globalCells, mesh_, zoneID, allStencil, allCellTypes);
 
    if (debug)
    {
        tmp<volScalarField> tfld
        (
            createField(mesh_, "allCellTypes_hole", allCellTypes)
        );
        tfld().write();
    }
    if (debug)
    {
        labelList stencilSize(mesh_.nCells());
        forAll(allStencil, celli)
        {
            stencilSize[celli] = allStencil[celli].size();
        }
        tmp<volScalarField> tfld
        (
            createField(mesh_, "allStencil_hole", stencilSize)
        );
        tfld().write();
    }
 
 
    // Add buffer interpolation layer(s) around holes
    scalarField allWeight(mesh_.nCells(), Zero);
    walkFront(layerRelax, allStencil, allCellTypes, allWeight);
 
    if (debug)
    {
        tmp<volScalarField> tfld
        (
            createField(mesh_, "allCellTypes_front", allCellTypes)
        );
        tfld().write();
    }
 
 
    // Convert cell-cell addressing to stencil in compact notation
 
    cellTypes_.transfer(allCellTypes);
    cellStencil_.setSize(mesh_.nCells());
    cellInterpolationWeights_.setSize(mesh_.nCells());
    DynamicList<label> interpolationCells;
    forAll(cellTypes_, cellI)
    {
        if (cellTypes_[cellI] == INTERPOLATED)
        {
            cellStencil_[cellI].transfer(allStencil[cellI]);
            cellInterpolationWeights_[cellI].setSize(1);
            cellInterpolationWeights_[cellI][0] = 1.0;
            interpolationCells.append(cellI);
        }
        else
        {
            cellStencil_[cellI].clear();
            cellInterpolationWeights_[cellI].clear();
        }
    }
    interpolationCells_.transfer(interpolationCells);
 
    List<Map<label>> compactMap;
    cellInterpolationMap_.reset
    (
        new mapDistribute(globalCells, cellStencil_, compactMap)
    );
    cellInterpolationWeight_.transfer(allWeight);
    dynamicOversetFvMesh::correctBoundaryConditions
    <
        volScalarField,
        oversetFvPatchField<scalar>
    >(cellInterpolationWeight_.boundaryFieldRef(), false);
 
 
    if (debug&2)
    {
        // Dump mesh
        mesh_.time().write();
 
        // Dump stencil
        mkDir(mesh_.time().timePath());
        OBJstream str(mesh_.time().timePath()/"injectionStencil.obj");
        Pout<< type() << " : dumping injectionStencil to "
            << str.name() << endl;
        pointField cc(mesh_.cellCentres());
        cellInterpolationMap().distribute(cc);
 
        forAll(cellStencil_, celli)
        {
            const labelList& slots = cellStencil_[celli];
            if (slots.size())
            {
                const point& accCc = mesh_.cellCentres()[celli];
                forAll(slots, i)
                {
                    const point& donorCc = cc[slots[i]];
                    str.write(linePointRef(accCc, 0.1*accCc+0.9*donorCc));
                }
            }
        }
    }
 
 
    // Extend stencil to get inverse distance weighted neighbours
    createStencil(globalCells);
 
 
    if (debug&2)
    {
        // Dump weight
        cellInterpolationWeight_.instance() = mesh_.time().timeName();
        cellInterpolationWeight_.write();
 
        // Dump max weight
        {
            scalarField maxMagWeight(mesh_.nCells(), Zero);
            forAll(cellStencil_, celli)
            {
                const scalarList& wghts = cellInterpolationWeights_[celli];
                forAll(wghts, i)
                {
                    if (mag(wghts[i]) > mag(maxMagWeight[celli]))
                    {
                        maxMagWeight[celli] = wghts[i];
                    }
                }
                if (mag(maxMagWeight[celli]) > 1)
                {
                    const pointField& cc = mesh_.cellCentres();
                    Pout<< "cell:" << celli
                        << " at:" << cc[celli]
                        << " zone:" << zoneID[celli]
                        << " donors:" << cellStencil_[celli]
                        << " weights:" << wghts
                        << " coords:"
                        << UIndirectList<point>(cc, cellStencil_[celli])
                        << " donorZone:"
                        << UIndirectList<label>(zoneID, cellStencil_[celli])
                        << endl;
                }
            }
            tmp<volScalarField> tfld
            (
                createField(mesh_, "maxMagWeight", maxMagWeight)
            );
            dynamicOversetFvMesh::correctBoundaryConditions
            <
                volScalarField,
                oversetFvPatchField<scalar>
            >(tfld.ref().boundaryFieldRef(), false);
            tfld().write();
        }
 
        // Dump cell types
        {
            tmp<volScalarField> tfld
            (
                createField(mesh_, "cellTypes", cellTypes_)
            );
            //tfld.ref().correctBoundaryConditions();
            dynamicOversetFvMesh::correctBoundaryConditions
            <
                volScalarField,
                oversetFvPatchField<scalar>
            >(tfld.ref().boundaryFieldRef(), false);
            tfld().write();
        }
 
 
        // Dump stencil
        mkDir(mesh_.time().timePath());
        OBJstream str(mesh_.time().timePath()/"stencil.obj");
        Pout<< type() << " : dumping to " << str.name() << endl;
        pointField cc(mesh_.cellCentres());
        cellInterpolationMap().distribute(cc);
 
        forAll(cellStencil_, celli)
        {
            const labelList& slots = cellStencil_[celli];
            if (slots.size())
            {
                const point& accCc = mesh_.cellCentres()[celli];
                forAll(slots, i)
                {
                    const point& donorCc = cc[slots[i]];
                    str.write(linePointRef(accCc, 0.1*accCc+0.9*donorCc));
                }
            }
        }
    }
 
    // Print some stats
    {
        labelList nCells(count(3, cellTypes_));
 
        label nLocal = 0;
        label nMixed = 0;
        label nRemote = 0;
        forAll(interpolationCells_, i)
        {
            label celli = interpolationCells_[i];
            const labelList& slots = cellStencil_[celli];
 
            bool hasLocal = false;
            bool hasRemote = false;
 
            forAll(slots, sloti)
            {
                if (slots[sloti] >= mesh_.nCells())
                {
                    hasRemote = true;
                }
                else
                {
                    hasLocal = true;
                }
            }
 
            if (hasRemote)
            {
                if (!hasLocal)
                {
                    nRemote++;
                }
                else
                {
                    nMixed++;
                }
            }
            else if (hasLocal)
            {
                nLocal++;
            }
        }
        reduce(nLocal, sumOp<label>());
        reduce(nMixed, sumOp<label>());
        reduce(nRemote, sumOp<label>());
 
        Info<< "Overset analysis : nCells : "
            << returnReduce(cellTypes_.size(), sumOp<label>()) << nl
            << incrIndent
            << indent << "calculated   : " << nCells[CALCULATED] << nl
            << indent << "interpolated : " << nCells[INTERPOLATED]
            << " (interpolated from local:" << nLocal
            << "  mixed local/remote:" << nMixed
            << "  remote:" << nRemote << ")" << nl
            << indent << "hole         : " << nCells[HOLE] << nl
            << decrIndent << endl;
    }
 
    // Tbd: detect if anything changed. Most likely it did!
    return true;
}
 
 
// ************************************************************************* //