From 1c64911367aea41033b98bcc99b22d1cf6013885 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Wed, 20 Jun 2018 07:57:25 +0100
Subject: [PATCH] ENH: Updated meshToMesh to be able to use the new
 processorLOD

---
 .../preProcessing/mapFieldsPar/mapFieldsPar.C |  48 ++-
 .../interRegionOption/interRegionOption.C     |  13 +-
 .../cellVolumeWeightCellCellStencil.C         |   3 +-
 .../meshToMeshMethod/meshToMeshMethod.C       |  37 ++-
 src/sampling/meshToMesh/meshToMesh.C          |  29 +-
 src/sampling/meshToMesh/meshToMesh.H          |  20 +-
 .../meshToMesh/meshToMeshParallelOps.C        | 277 ++++++++++--------
 7 files changed, 281 insertions(+), 146 deletions(-)

diff --git a/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
index 1d021b5280..a44bcb651c 100644
--- a/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
+++ b/applications/utilities/preProcessing/mapFieldsPar/mapFieldsPar.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,6 +46,7 @@ void mapConsistentMesh
     const fvMesh& meshTarget,
     const word& mapMethod,
     const word& AMIMapMethod,
+    const word& procMapMethod,
     const bool subtract,
     const wordHashSet& selectedFields,
     const bool noLagrangian
@@ -54,7 +55,14 @@ void mapConsistentMesh
     Info<< nl << "Consistently creating and mapping fields for time "
         << meshSource.time().timeName() << nl << endl;
 
-    meshToMesh interp(meshSource, meshTarget, mapMethod, AMIMapMethod);
+    meshToMesh interp
+    (
+        meshSource,
+        meshTarget,
+        mapMethod,
+        AMIMapMethod,
+        meshToMesh::procMapMethodNames_[procMapMethod]
+    );
 
     if (subtract)
     {
@@ -85,6 +93,7 @@ void mapSubMesh
     const wordList& cuttingPatches,
     const word& mapMethod,
     const word& AMIMapMethod,
+    const word& procMapMethod,
     const bool subtract,
     const wordHashSet& selectedFields,
     const bool noLagrangian
@@ -100,7 +109,8 @@ void mapSubMesh
         mapMethod,
         AMIMapMethod,
         patchMap,
-        cuttingPatches
+        cuttingPatches,
+        meshToMesh::procMapMethodNames_[procMapMethod]
     );
 
     if (subtract)
@@ -171,6 +181,12 @@ int main(int argc, char *argv[])
         "word",
         "specify the patch mapping method (direct|mapNearest|faceAreaWeight)"
     );
+    argList::addOption
+    (
+        "procMapMethod",
+        "word",
+        "specify the processor distribution map method (AABB|LOD)"
+    );
     argList::addBoolOption
     (
         "subtract",
@@ -217,7 +233,7 @@ int main(int argc, char *argv[])
 
     word mapMethod = meshToMesh::interpolationMethodNames_
     [
-        meshToMesh::imCellVolumeWeight
+        meshToMesh::interpolationMethod::imCellVolumeWeight
     ];
 
     if  (args.readIfPresent("mapMethod", mapMethod))
@@ -225,7 +241,6 @@ int main(int argc, char *argv[])
         Info<< "Mapping method: " << mapMethod << endl;
     }
 
-
     word patchMapMethod;
     if (meshToMesh::interpolationMethodNames_.found(mapMethod))
     {
@@ -233,12 +248,25 @@ int main(int argc, char *argv[])
         meshToMesh::interpolationMethod method =
             meshToMesh::interpolationMethodNames_[mapMethod];
 
-        patchMapMethod = AMIPatchToPatchInterpolation::interpolationMethodNames_
-        [
-            meshToMesh::interpolationMethodAMI(method)
-        ];
+        patchMapMethod =
+            AMIPatchToPatchInterpolation::interpolationMethodNames_
+            [
+                meshToMesh::interpolationMethodAMI(method)
+            ];
     }
 
+    word procMapMethod =
+        meshToMesh::procMapMethodNames_
+        [
+            meshToMesh::procMapMethod::pmAABB
+        ];
+
+    if (args.readIfPresent("prcoMapMethod", procMapMethod))
+    {
+        Info<< "Processor map method: " << procMapMethod << endl;
+    }
+
+
     // Optionally override
     if (args.readIfPresent("patchMapMethod", patchMapMethod))
     {
@@ -325,6 +353,7 @@ int main(int argc, char *argv[])
             meshTarget,
             mapMethod,
             patchMapMethod,
+            procMapMethod,
             subtract,
             selectedFields,
             noLagrangian
@@ -340,6 +369,7 @@ int main(int argc, char *argv[])
             cuttingPatches,
             mapMethod,
             patchMapMethod,
+            procMapMethod,
             subtract,
             selectedFields,
             noLagrangian
diff --git a/src/fvOptions/interRegionOption/interRegionOption.C b/src/fvOptions/interRegionOption/interRegionOption.C
index 32a3f88564..bf81ccf383 100644
--- a/src/fvOptions/interRegionOption/interRegionOption.C
+++ b/src/fvOptions/interRegionOption/interRegionOption.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,10 +65,17 @@ void Foam::fv::interRegionOption::setMapper()
                 (
                     mesh_,
                     nbrMesh,
-                    meshToMesh::interpolationMethodNames_.lookup
+                    meshToMesh::interpolationMethodNames_.lookupOrDefault
                     (
                         "interpolationMethod",
-                        coeffs_
+                        coeffs_,
+                        meshToMesh::interpolationMethod::imCellVolumeWeight
+                    ),
+                    meshToMesh::procMapMethodNames_.lookupOrDefault
+                    (
+                        "procMapMethod",
+                        coeffs_,
+                        meshToMesh::procMapMethod::pmAABB
                     ),
                     false // not interpolating patches
                 )
diff --git a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
index 8ffc7f6c6e..1c9a5e72c4 100644
--- a/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
+++ b/src/overset/cellCellStencil/cellVolumeWeight/cellVolumeWeightCellCellStencil.C
@@ -784,9 +784,10 @@ bool Foam::cellCellStencils::cellVolumeWeight::update()
             (
                 srcMesh,
                 tgtMesh,
-                meshToMesh::imCellVolumeWeight,
+                meshToMesh::interpolationMethod::imCellVolumeWeight,
                 HashTable<word>(0),     // patchMap,
                 wordList(0),            // cuttingPatches
+                meshToMesh::procMapMethod::pmAABB,
                 false                   // do not normalise
             );
 
diff --git a/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
index afb6b727dc..ef21a395f1 100644
--- a/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
+++ b/src/sampling/meshToMesh/calcMethod/meshToMeshMethod/meshToMeshMethod.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
-     \\/     M anipulation  |
+     \\/     M anipulation  | Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -84,7 +84,14 @@ bool Foam::meshToMeshMethod::intersect
 
     tetOverlapVolume overlapEngine;
 
-    treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
+    // Note: avoid demand-driven construction of cellPoints
+    // treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
+    const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
+    treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
+    for (label i = 1; i < cellFaces.size(); ++i)
+    {
+        bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
+    }
 
     return overlapEngine.cellCellOverlapMinDecomp
     (
@@ -106,7 +113,14 @@ Foam::scalar Foam::meshToMeshMethod::interVol
 {
     tetOverlapVolume overlapEngine;
 
-    treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
+    // Note: avoid demand-driven construction of cellPoints
+    // treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
+    const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
+    treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
+    for (label i = 1; i < cellFaces.size(); ++i)
+    {
+        bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
+    }
 
     scalar vol = overlapEngine.cellCellOverlapVolumeMinDecomp
     (
@@ -124,21 +138,28 @@ Foam::scalar Foam::meshToMeshMethod::interVol
 Foam::Tuple2<Foam::scalar, Foam::point>
 Foam::meshToMeshMethod::interVolAndCentroid
 (
-    const label srcCellI,
-    const label tgtCellI
+    const label srcCelli,
+    const label tgtCelli
 )
 {
     tetOverlapVolume overlapEngine;
 
-    treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCellI]);
+    // Note: avoid demand-driven construction of cellPoints
+    // treeBoundBox bbTgtCell(tgt_.points(), tgt_.cellPoints()[tgtCelli]);
+    const UList<label>& cellFaces = tgt_.cells()[tgtCelli];
+    treeBoundBox bbTgtCell(tgt_.points(), tgt_.faces()[cellFaces[0]]);
+    for (label i = 1; i < cellFaces.size(); ++i)
+    {
+        bbTgtCell.add(tgt_.points(), tgt_.faces()[cellFaces[i]]);
+    }
 
     Tuple2<scalar, point> volAndInertia =
     overlapEngine.cellCellOverlapMomentMinDecomp
     (
         src_,
-        srcCellI,
+        srcCelli,
         tgt_,
-        tgtCellI,
+        tgtCelli,
         bbTgtCell
     );
 
diff --git a/src/sampling/meshToMesh/meshToMesh.C b/src/sampling/meshToMesh/meshToMesh.C
index d4ab2b7449..2c6a5e8ab6 100644
--- a/src/sampling/meshToMesh/meshToMesh.C
+++ b/src/sampling/meshToMesh/meshToMesh.C
@@ -52,6 +52,17 @@ Foam::meshToMesh::interpolationMethodNames_
 };
 
 
+const Foam::Enum
+<
+    Foam::meshToMesh::procMapMethod
+>
+Foam::meshToMesh::procMapMethodNames_
+{
+    { procMapMethod::pmAABB, "AABB" },
+    { procMapMethod::pmLOD, "LOD" },
+};
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<>
@@ -628,18 +639,18 @@ Foam::meshToMesh::interpolationMethodAMI(const interpolationMethod method)
 {
     switch (method)
     {
-        case imDirect:
+        case interpolationMethod::imDirect:
         {
             return AMIPatchToPatchInterpolation::imDirect;
             break;
         }
-        case imMapNearest:
+        case interpolationMethod::imMapNearest:
         {
             return AMIPatchToPatchInterpolation::imMapNearest;
             break;
         }
-        case imCellVolumeWeight:
-        case imCorrectedCellVolumeWeight:
+        case interpolationMethod::imCellVolumeWeight:
+        case interpolationMethod::imCorrectedCellVolumeWeight:
         {
             return AMIPatchToPatchInterpolation::imFaceAreaWeight;
             break;
@@ -647,7 +658,7 @@ Foam::meshToMesh::interpolationMethodAMI(const interpolationMethod method)
         default:
         {
             FatalErrorInFunction
-                << "Unhandled enumeration " << method
+                << "Unhandled enumeration " << interpolationMethodNames_[method]
                 << abort(FatalError);
         }
     }
@@ -827,11 +838,13 @@ Foam::meshToMesh::meshToMesh
     const polyMesh& src,
     const polyMesh& tgt,
     const interpolationMethod& method,
+    const procMapMethod& mapMethod,
     bool interpAllPatches
 )
 :
     srcRegion_(src),
     tgtRegion_(tgt),
+    procMapMethod_(mapMethod),
     srcPatchID_(),
     tgtPatchID_(),
     patchAMIs_(),
@@ -865,11 +878,13 @@ Foam::meshToMesh::meshToMesh
     const polyMesh& tgt,
     const word& methodName,
     const word& AMIMethodName,
+    const procMapMethod& mapMethod,
     bool interpAllPatches
 )
 :
     srcRegion_(src),
     tgtRegion_(tgt),
+    procMapMethod_(mapMethod),
     srcPatchID_(),
     tgtPatchID_(),
     patchAMIs_(),
@@ -896,11 +911,13 @@ Foam::meshToMesh::meshToMesh
     const interpolationMethod& method,
     const HashTable<word>& patchMap,
     const wordList& cuttingPatches,
+    const procMapMethod& mapMethod,
     const bool normalise
 )
 :
     srcRegion_(src),
     tgtRegion_(tgt),
+    procMapMethod_(mapMethod),
     srcPatchID_(),
     tgtPatchID_(),
     patchAMIs_(),
@@ -936,11 +953,13 @@ Foam::meshToMesh::meshToMesh
     const word& AMIMethodName,  // boundary mapping
     const HashTable<word>& patchMap,
     const wordList& cuttingPatches,
+    const procMapMethod& mapMethod,
     const bool normalise
 )
 :
     srcRegion_(src),
     tgtRegion_(tgt),
+    procMapMethod_(mapMethod),
     srcPatchID_(),
     tgtPatchID_(),
     patchAMIs_(),
diff --git a/src/sampling/meshToMesh/meshToMesh.H b/src/sampling/meshToMesh/meshToMesh.H
index e1567c18a4..40aaaf329a 100644
--- a/src/sampling/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMesh/meshToMesh.H
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2012-2016 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2016 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,7 +66,7 @@ public:
     // Public data types
 
         //- Enumeration specifying interpolation method
-        enum interpolationMethod
+        enum class interpolationMethod
         {
             imDirect,
             imMapNearest,
@@ -76,6 +76,15 @@ public:
 
         static const Enum<interpolationMethod> interpolationMethodNames_;
 
+        //- Enumeration specifying processor parallel map construction method
+        enum class procMapMethod
+        {
+            pmAABB,
+            pmLOD
+        };
+
+        static const Enum<procMapMethod> procMapMethodNames_;
+
 private:
 
     // Private data
@@ -86,6 +95,9 @@ private:
         //- Reference to the target mesh
         const polyMesh& tgtRegion_;
 
+        //- Processor parallel map construction method
+        procMapMethod procMapMethod_;
+
         //- List of target patch IDs per source patch (local index)
         List<label> srcPatchID_;
 
@@ -301,6 +313,7 @@ public:
         const polyMesh& src,
         const polyMesh& tgt,
         const interpolationMethod& method,
+        const procMapMethod& mapMethod = procMapMethod::pmAABB,
         const bool interpAllPatches = true
     );
 
@@ -311,6 +324,7 @@ public:
         const polyMesh& tgt,
         const word& methodName,     // internal mapping
         const word& AMIMethodName,  // boundary mapping
+        const procMapMethod& mapMethod = procMapMethod::pmAABB,
         const bool interpAllPatches = true
     );
 
@@ -322,6 +336,7 @@ public:
         const interpolationMethod& method,
         const HashTable<word>& patchMap,
         const wordList& cuttingPatches,
+        const procMapMethod& mapMethod = procMapMethod::pmAABB,
         const bool normalise = true
     );
 
@@ -335,6 +350,7 @@ public:
         const word& AMIMethodName,  // boundary mapping
         const HashTable<word>& patchMap,
         const wordList& cuttingPatches,
+        const procMapMethod& mapMethod = procMapMethod::pmAABB,
         const bool normalise = true
     );
 
diff --git a/src/sampling/meshToMesh/meshToMeshParallelOps.C b/src/sampling/meshToMesh/meshToMeshParallelOps.C
index bae3aa442a..b391a6132f 100644
--- a/src/sampling/meshToMesh/meshToMeshParallelOps.C
+++ b/src/sampling/meshToMesh/meshToMeshParallelOps.C
@@ -3,7 +3,7 @@
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
     \\  /    A nd           | Copyright (C) 2012-2017 OpenFOAM Foundation
-     \\/     M anipulation  | Copyright (C) 2015-2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2015-2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,6 +31,7 @@ License
 #include "processorPolyPatch.H"
 #include "SubField.H"
 #include "AABBTree.H"
+#include "cellBox.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -98,12 +99,12 @@ Foam::label Foam::meshToMesh::calcOverlappingProcs
     {
         const treeBoundBoxList& bbp = procBb[proci];
 
-        forAll(bbp, bbi)
+        for (const treeBoundBox& b : bbp)
         {
-            if (bbp[bbi].overlaps(bb))
+            if (b.overlaps(bb))
             {
                 overlaps[proci] = true;
-                nOverlaps++;
+                ++nOverlaps;
                 break;
             }
         }
@@ -119,139 +120,179 @@ Foam::autoPtr<Foam::mapDistribute> Foam::meshToMesh::calcProcMap
     const polyMesh& tgt
 ) const
 {
-    // get decomposition of cells on src mesh
-    List<treeBoundBoxList> procBb(Pstream::nProcs());
-
-    if (src.nCells() > 0)
+    switch (procMapMethod_)
     {
-        procBb[Pstream::myProcNo()] = AABBTree<labelList>
-        (
-            src.cellPoints(),
-            src.points(),
-            false
-        ).boundBoxes();
-    }
-    else
-    {
-        procBb[Pstream::myProcNo()] = treeBoundBoxList();
-    }
-
-
-    Pstream::gatherList(procBb);
-    Pstream::scatterList(procBb);
-
-
-    if (debug)
-    {
-        InfoInFunction
-            << "Determining extent of src mesh per processor:" << nl
-            << "\tproc\tbb" << endl;
-        forAll(procBb, proci)
+        case procMapMethod::pmLOD:
         {
-            Info<< '\t' << proci << '\t' << procBb[proci] << endl;
+            Info<< "meshToMesh: Using processorLOD method" << endl;
+
+            // Create processor map of overlapping faces. This map gets
+            // (possibly remote) cells from the tgt mesh such that they
+            // (together) cover all of the src mesh
+            label nGlobalSrcCells =
+                returnReduce(src.cells().size(), sumOp<label>());
+            label cellsPerBox = max(1, 0.001*nGlobalSrcCells);
+            typename processorLODs::cellBox boxLOD
+            (
+                src.cells(),
+                src.faces(),
+                src.points(),
+                tgt.cells(),
+                tgt.faces(),
+                tgt.points(),
+                cellsPerBox,
+                src.nCells()
+            );
+
+            return boxLOD.map();
+            break;
         }
-    }
-
-
-    // determine which cells of tgt mesh overlaps src mesh per proc
-    const cellList& cells = tgt.cells();
-    const faceList& faces = tgt.faces();
-    const pointField& points = tgt.points();
-
-    labelListList sendMap;
-
-    {
-        // per processor indices into all segments to send
-        List<DynamicList<label>> dynSendMap(Pstream::nProcs());
-        label iniSize = floor(tgt.nCells()/Pstream::nProcs());
-
-        forAll(dynSendMap, proci)
+        default:
         {
-            dynSendMap[proci].setCapacity(iniSize);
-        }
+            Info<< "meshToMesh: Using AABBTree method" << endl;
 
-        // work array - whether src processor bb overlaps the tgt cell bounds
-        boolList procBbOverlaps(Pstream::nProcs());
-        forAll(cells, celli)
-        {
-            const cell& c = cells[celli];
+            // get decomposition of cells on src mesh
+            List<treeBoundBoxList> procBb(Pstream::nProcs());
 
-            // determine bounding box of tgt cell
-            boundBox cellBb(boundBox::invertedBox);
-            forAll(c, facei)
+            if (src.nCells() > 0)
             {
-                const face& f = faces[c[facei]];
-                forAll(f, fp)
+                procBb[Pstream::myProcNo()] = AABBTree<labelList>
+                (
+                    src.cellPoints(),
+                    src.points(),
+                    false
+                ).boundBoxes();
+            }
+            else
+            {
+                procBb[Pstream::myProcNo()] = treeBoundBoxList();
+            }
+
+
+            Pstream::gatherList(procBb);
+            Pstream::scatterList(procBb);
+
+
+            if (debug)
+            {
+                InfoInFunction
+                    << "Determining extent of src mesh per processor:" << nl
+                    << "\tproc\tbb" << endl;
+                forAll(procBb, proci)
                 {
-                    cellBb.add(points, f);
+                    Info<< '\t' << proci << '\t' << procBb[proci] << endl;
                 }
             }
 
-            // find the overlapping tgt cells on each src processor
-            (void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);
 
-            forAll(procBbOverlaps, proci)
+            // determine which cells of tgt mesh overlaps src mesh per proc
+            const cellList& cells = tgt.cells();
+            const faceList& faces = tgt.faces();
+            const pointField& points = tgt.points();
+
+            labelListList sendMap;
+
             {
-                if (procBbOverlaps[proci])
+                // per processor indices into all segments to send
+                List<DynamicList<label>> dynSendMap(Pstream::nProcs());
+                label iniSize = floor(tgt.nCells()/Pstream::nProcs());
+
+                forAll(dynSendMap, proci)
                 {
-                    dynSendMap[proci].append(celli);
+                    dynSendMap[proci].setCapacity(iniSize);
+                }
+
+                // work array - whether src processor bb overlaps the tgt cell
+                // bounds
+                boolList procBbOverlaps(Pstream::nProcs());
+                forAll(cells, celli)
+                {
+                    const cell& c = cells[celli];
+
+                    // determine bounding box of tgt cell
+                    boundBox cellBb(boundBox::invertedBox);
+                    forAll(c, facei)
+                    {
+                        const face& f = faces[c[facei]];
+                        forAll(f, fp)
+                        {
+                            cellBb.add(points, f);
+                        }
+                    }
+
+                    // find the overlapping tgt cells on each src processor
+                    (void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);
+
+                    forAll(procBbOverlaps, proci)
+                    {
+                        if (procBbOverlaps[proci])
+                        {
+                            dynSendMap[proci].append(celli);
+                        }
+                    }
+                }
+
+                // convert dynamicList to labelList
+                sendMap.setSize(Pstream::nProcs());
+                forAll(sendMap, proci)
+                {
+                    sendMap[proci].transfer(dynSendMap[proci]);
                 }
             }
-        }
 
-        // convert dynamicList to labelList
-        sendMap.setSize(Pstream::nProcs());
-        forAll(sendMap, proci)
-        {
-            sendMap[proci].transfer(dynSendMap[proci]);
+            // debug printing
+            if (debug)
+            {
+                Pout<< "Of my " << cells.size()
+                    << " target cells I need to send to:" << nl
+                    << "\tproc\tcells" << endl;
+                forAll(sendMap, proci)
+                {
+                    Pout<< '\t' << proci << '\t' << sendMap[proci].size()
+                        << endl;
+                }
+            }
+
+
+            // send over how many tgt cells I need to receive from each
+            // processor
+            labelListList sendSizes(Pstream::nProcs());
+            sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
+            forAll(sendMap, proci)
+            {
+                sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
+            }
+            Pstream::gatherList(sendSizes);
+            Pstream::scatterList(sendSizes);
+
+
+            // determine order of receiving
+            labelListList constructMap(Pstream::nProcs());
+
+            label segmentI = 0;
+            forAll(constructMap, proci)
+            {
+                // what I need to receive is what other processor is sending
+                // to me
+                label nRecv = sendSizes[proci][Pstream::myProcNo()];
+                constructMap[proci].setSize(nRecv);
+
+                for (label i = 0; i < nRecv; i++)
+                {
+                    constructMap[proci][i] = segmentI++;
+                }
+            }
+
+            return autoPtr<mapDistribute>::New
+            (
+                segmentI,       // size after construction
+                std::move(sendMap),
+                std::move(constructMap)
+            );
+
+            break;
         }
     }
-
-    // debug printing
-    if (debug)
-    {
-        Pout<< "Of my " << cells.size() << " target cells I need to send to:"
-            << nl << "\tproc\tcells" << endl;
-        forAll(sendMap, proci)
-        {
-            Pout<< '\t' << proci << '\t' << sendMap[proci].size() << endl;
-        }
-    }
-
-
-    // send over how many tgt cells I need to receive from each processor
-    labelListList sendSizes(Pstream::nProcs());
-    sendSizes[Pstream::myProcNo()].setSize(Pstream::nProcs());
-    forAll(sendMap, proci)
-    {
-        sendSizes[Pstream::myProcNo()][proci] = sendMap[proci].size();
-    }
-    Pstream::gatherList(sendSizes);
-    Pstream::scatterList(sendSizes);
-
-
-    // determine order of receiving
-    labelListList constructMap(Pstream::nProcs());
-
-    label segmentI = 0;
-    forAll(constructMap, proci)
-    {
-        // what I need to receive is what other processor is sending to me
-        label nRecv = sendSizes[proci][Pstream::myProcNo()];
-        constructMap[proci].setSize(nRecv);
-
-        for (label i = 0; i < nRecv; i++)
-        {
-            constructMap[proci][i] = segmentI++;
-        }
-    }
-
-    return autoPtr<mapDistribute>::New
-    (
-        segmentI,       // size after construction
-        std::move(sendMap),
-        std::move(constructMap)
-    );
 }