From 33527b9bd3813a788748bda4237d1bbbb9aa4741 Mon Sep 17 00:00:00 2001
From: Sopheak Seng <sopheak.seng@bureauveritas.com>
Date: Tue, 18 Jul 2017 18:56:12 +0200
Subject: [PATCH 1/3] BUGFIX: when run in parallel, refineMesh needs "cellCuts"
 to sync. processor patches

---
 src/dynamicMesh/meshCut/cellCuts/cellCuts.C     | 141 ++++++++++++++++++++++++
 src/dynamicMesh/meshCut/cellCuts/cellCuts.H     |   2 +
 src/dynamicMesh/meshCut/directions/directions.C |   2 +-
 3 files changed, 144 insertions(+), 1 deletion(-)

diff --git a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
index b49a709..efca831 100644
--- a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
+++ b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
@@ -34,6 +34,9 @@ License
 #include "OFstream.H"
 #include "plane.H"
 
+#include "processorPolyPatch.H"
+#include "EdgeMap.H"
+
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
@@ -114,6 +117,142 @@ Foam::label Foam::cellCuts::firstUnique
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::cellCuts::syncProc(const polyMesh& mesh)
+{
+    if (!Pstream::parRun()) return;
+
+    forAll(pointIsCut_, i)
+    {
+        if (pointIsCut_[i])
+        {
+            WarningInFunction
+            << "CutPoints parallel sync. is not yet implemented ... "
+            << endl;
+            break;
+        }
+    }
+
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+    const polyBoundaryMesh& bm(mesh.boundaryMesh());
+    const labelListList& faceEdges(mesh.faceEdges());
+
+    // Send edgeIsCut_
+    forAll(bm, i)
+    {
+        if (!isA<processorPolyPatch>(bm[i])) continue;
+        const processorPolyPatch& pp =
+                refCast<const processorPolyPatch>(bm[i]);
+        const edgeList& edges(pp.edges());
+        const labelList& meshEdges(pp.meshEdges());
+        const labelList& nbrPts(pp.neighbPoints());
+        EdgeMap<scalar> cutEdges(meshEdges.size());
+        forAll(meshEdges, j)
+        {
+            label ei(meshEdges[j]);
+            if (!edgeIsCut_[ei]) continue;
+            const edge& e(edges[j]);
+            const edge nbrEdge(nbrPts[e[0]], nbrPts[e[1]]);
+            cutEdges.insert(nbrEdge, edgeWeight_[ei]);
+        }
+        UOPstream toNeighb(pp.neighbProcNo(), pBufs);
+        toNeighb << cutEdges;
+    }
+    pBufs.finishedSends();
+    // Receive edgeIsCut_ from neighbour(s)
+    forAll(bm, i)
+    {
+        if (!isA<processorPolyPatch>(bm[i])) continue;
+        const processorPolyPatch& pp =
+                refCast<const processorPolyPatch>(bm[i]);
+        const labelList& meshEdges(pp.meshEdges());
+        EdgeMap<scalar> nbrCutEdges;
+        {
+            UIPstream fromNbr(pp.neighbProcNo(), pBufs);
+            fromNbr >> nbrCutEdges;
+        }
+        forAllConstIter(EdgeMap<scalar>, nbrCutEdges, iter)
+        {
+            label ei(meshEdges[pp.whichEdge(iter.key())]);
+            if (edgeIsCut_[ei])
+            {
+                if (fabs(edgeWeight_[ei] - iter()) > SMALL)
+                {
+                    FatalErrorInFunction
+                    << "Same edges are cut differently on proc(s) ..."
+                    << abort(FatalError);
+                }
+            }
+            else
+            {
+                edgeIsCut_[ei] = true;
+                edgeWeight_[ei] = iter();
+            }
+        }
+    }
+
+    // count new point(s) on edgeIsCut_
+    labelList addedPointLabels(edgeIsCut_.size(),-1);
+    forAll(edgeIsCut_, i)
+    {
+        if (edgeIsCut_[i]) addedPointLabels[i] = mesh.nPoints() + i;
+    }
+
+    // Send faceSplitCut_
+    forAll(bm, i)
+    {
+        if (!isA<processorPolyPatch>(bm[i])) continue;
+        const processorPolyPatch& pp =
+                refCast<const processorPolyPatch>(bm[i]);
+        label start(pp.start());
+        Map<edge> myFaceSplit(pp.size());
+        forAll(pp, k)
+        {
+            label fi(start + k);
+            Map<edge>::const_iterator iter(faceSplitCut_.find(fi));
+            if (iter == faceSplitCut_.end()) continue;
+            const edge& splitEdge(iter());
+
+            // FIXME: splitEdge not used on receive
+            myFaceSplit.insert(k, splitEdge);
+        }
+        UOPstream toNeighb(pp.neighbProcNo(), pBufs);
+        toNeighb << myFaceSplit;
+    }
+    pBufs.finishedSends();
+    // Receive faceSplitCut_ from neighbour(s)
+    forAll(bm, i)
+    {
+        if (!isA<processorPolyPatch>(bm[i])) continue;
+        const processorPolyPatch& pp =
+                refCast<const processorPolyPatch>(bm[i]);
+        label start(pp.start());
+        Map<edge> nbrFaceSplit;
+        {
+            UIPstream fromNbr(pp.neighbProcNo(), pBufs);
+            fromNbr >> nbrFaceSplit;
+        }
+        forAllConstIter(Map<edge>, nbrFaceSplit, iter)
+        {
+            label fi(start + iter.key());
+            bool found(faceSplitCut_.find(fi) != faceSplitCut_.end());
+            if (found) continue;
+            const labelList& fe(faceEdges[fi]);
+            DynamicList<label> cuts(fe.size());
+            forAll(fe, k)
+            {
+                if (edgeIsCut_[fe[k]]) cuts.append(addedPointLabels[fe[k]]);
+            }
+            if (cuts.size()!=2)
+            {
+                FatalErrorInFunction
+                << "Face split needs 2 edges ... " << abort(FatalError);
+            }
+            const edge e(cuts[0], cuts[1]);
+            faceSplitCut_.insert(fi, e);
+        }
+    }
+}
+
 void Foam::cellCuts::writeUncutOBJ
 (
     const fileName& dir,
@@ -2831,6 +2970,8 @@ Foam::cellCuts::cellCuts
 
     clearOut();
 
+    syncProc(mesh);
+
     if (debug)
     {
         Pout<< "cellCuts : leaving constructor from cellCutter" << endl;
diff --git a/src/dynamicMesh/meshCut/cellCuts/cellCuts.H b/src/dynamicMesh/meshCut/cellCuts/cellCuts.H
index 7e80ce1..c9da2f8 100644
--- a/src/dynamicMesh/meshCut/cellCuts/cellCuts.H
+++ b/src/dynamicMesh/meshCut/cellCuts/cellCuts.H
@@ -445,6 +445,8 @@ class cellCuts
             //- Check various consistencies.
             void check() const;
 
+            //- Sync. processor patches
+            void syncProc(const polyMesh& mesh);
 
         //- Disallow default bitwise copy construct
         cellCuts(const cellCuts&);
diff --git a/src/dynamicMesh/meshCut/directions/directions.C b/src/dynamicMesh/meshCut/directions/directions.C
index 2097a2c..9d7f932 100644
--- a/src/dynamicMesh/meshCut/directions/directions.C
+++ b/src/dynamicMesh/meshCut/directions/directions.C
@@ -322,7 +322,7 @@ Foam::directions::directions
         vector normal = tan1 ^ tan2;
         normal /= mag(normal);
 
-        Pout<< "Global Coordinate system:" << endl
+        Info<< "Global Coordinate system:" << endl
             << "     normal : " << normal << endl
             << "     tan1   : " << tan1 << endl
             << "     tan2   : " << tan2
-- 
1.9.1


From c0c5b4d8d558d687950f845ae0f5704200796436 Mon Sep 17 00:00:00 2001
From: Sopheak Seng <sopheak.seng@bureauveritas.com>
Date: Wed, 19 Jul 2017 11:31:09 +0200
Subject: [PATCH 2/3] BUG: refineMesh: meshCutter should allow parallel run
 although there is no cell to cut

---
 src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C b/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C
index 903e1c6..8c22d2a 100644
--- a/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C
+++ b/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C
@@ -536,7 +536,7 @@ void Foam::meshCutter::setRefinement
     addedPoints_.clear();
     addedPoints_.resize(cuts.nLoops());
 
-    if (cuts.nLoops() == 0)
+    if ((cuts.nLoops() == 0) && (!Pstream::parRun()))
     {
         return;
     }
-- 
1.9.1


From f720bb21a5cc4368fc4d2394094a6f305a0b4088 Mon Sep 17 00:00:00 2001
From: Sopheak Seng <sopheak.seng@bureauveritas.com>
Date: Wed, 19 Jul 2017 11:32:32 +0200
Subject: [PATCH 3/3] BUG: refineMesh: edge statistic: should not count
 processor patches twice

---
 .../mesh/manipulation/refineMesh/refineMesh.C      | 105 ++++++++++++++-------
 1 file changed, 73 insertions(+), 32 deletions(-)

diff --git a/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C b/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
index debcd9f..7b33d90 100644
--- a/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
+++ b/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
@@ -52,6 +52,8 @@ Description
 #include "plane.H"
 #include "SubField.H"
 
+#include "processorPolyPatch.H"
+
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -61,72 +63,104 @@ static const scalar edgeTol = 1e-3;
 
 
 // Print edge statistics on mesh.
-void printEdgeStats(const primitiveMesh& mesh)
+//void printEdgeStats(const primitiveMesh& mesh)
+void printEdgeStats(const polyMesh& mesh)
 {
-    label nX = 0;
-    label nY = 0;
-    label nZ = 0;
+    label nX(0), nY(0), nZ(0), other(0);
 
-    scalar minX = GREAT;
-    scalar maxX = -GREAT;
+    scalar minX(GREAT), maxX(-GREAT);
+    scalar minY(GREAT), maxY(-GREAT);
+    scalar minZ(GREAT), maxZ(-GREAT);
+    scalar minOther(GREAT), maxOther(-GREAT);
     vector x(1, 0, 0);
-
-    scalar minY = GREAT;
-    scalar maxY = -GREAT;
     vector y(0, 1, 0);
-
-    scalar minZ = GREAT;
-    scalar maxZ = -GREAT;
     vector z(0, 0, 1);
 
-    scalar minOther = GREAT;
-    scalar maxOther = -GREAT;
-
     const edgeList& edges = mesh.edges();
 
+    // don't count processor patches twice(s)
+    label nXc(0), nYc(0), nZc(0), otherc(0);
+    const polyBoundaryMesh& bm(mesh.boundaryMesh());
+    const labelListList& edgeFaces(mesh.edgeFaces());
+
     forAll(edges, edgeI)
     {
-        const edge& e = edges[edgeI];
+        const labelList& faces = edgeFaces[edgeI];
+        bool onProcPatch(false);
+        forAll(faces, fi)
+        {
+            label patchi = bm.whichPatch(faces[fi]);
+            if (patchi<0) continue;
+            if (isA<processorPolyPatch>(bm[patchi]))
+            {
+                onProcPatch=true;
+                break;
+            }
+        }
 
+        const edge& e = edges[edgeI];
         vector eVec(e.vec(mesh.points()));
-
-        scalar eMag = mag(eVec);
-
-        eVec /= eMag;
-
+        scalar eMag = mag(eVec); eVec /= eMag;
         if (mag(eVec & x) > 1-edgeTol)
         {
             minX = min(minX, eMag);
             maxX = max(maxX, eMag);
             nX++;
+            if (onProcPatch) nXc++;
         }
         else if (mag(eVec & y) > 1-edgeTol)
         {
             minY = min(minY, eMag);
             maxY = max(maxY, eMag);
             nY++;
+            if (onProcPatch) nYc++;
         }
         else if (mag(eVec & z) > 1-edgeTol)
         {
             minZ = min(minZ, eMag);
             maxZ = max(maxZ, eMag);
             nZ++;
+            if (onProcPatch) nZc++;
         }
         else
         {
             minOther = min(minOther, eMag);
             maxOther = max(maxOther, eMag);
+            other++;
+            if (onProcPatch) otherc++;
         }
     }
 
-    Pout<< "Mesh edge statistics:" << endl
+    reduce(nX, sumOp<label>());
+    reduce(nY, sumOp<label>());
+    reduce(nZ, sumOp<label>());
+    reduce(other, sumOp<label>());
+    reduce(minX, minOp<scalar>());
+    reduce(minY, minOp<scalar>());
+    reduce(minZ, minOp<scalar>());
+    reduce(maxX, maxOp<scalar>());
+    reduce(maxY, maxOp<scalar>());
+    reduce(maxZ, maxOp<scalar>());
+    reduce(minOther, minOp<scalar>());
+    reduce(maxOther, maxOp<scalar>());
+    //
+    reduce(nXc, sumOp<label>());
+    reduce(nYc, sumOp<label>());
+    reduce(nZc, sumOp<label>());
+    reduce(otherc, sumOp<label>());
+    nX -= nXc/2;
+    nY -= nYc/2;
+    nZ -= nZc/2;
+    other -= otherc/2;
+
+    Info<< "Mesh edge statistics:" << endl
         << "    x aligned :  number:" << nX << "\tminLen:" << minX
         << "\tmaxLen:" << maxX << endl
         << "    y aligned :  number:" << nY << "\tminLen:" << minY
         << "\tmaxLen:" << maxY << endl
         << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
         << "\tmaxLen:" << maxZ << endl
-        << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
+        << "    other     :  number:" << other
         << "\tminLen:" << minOther
         << "\tmaxLen:" << maxOther << endl << endl;
 }
@@ -225,13 +259,20 @@ int main(int argc, char *argv[])
     {
         const word setName(refineDict.lookup("set"));
 
-        cellSet cells(mesh, setName);
-
-        Pout<< "Read " << cells.size() << " cells from cellSet "
-            << cells.instance()/cells.local()/cells.name()
-            << endl << endl;
+        label zoneId = mesh.cellZones().findZoneID(setName);
+        if (zoneId >= 0)
+        {
+            labelList cells(mesh.cellZones()[zoneId]);
+            refCells = cells;
+        }
+        else
+        {
+            cellSet cells(mesh, setName);
+            refCells = cells.toc();
+        }
 
-        refCells = cells.toc();
+        Info<< "Read " << returnReduce(refCells.size(), sumOp<scalar>())
+            << " cells from cellSet " << setName << endl << endl;
     }
     else
     {
@@ -314,13 +355,13 @@ int main(int argc, char *argv[])
 
 
     // Write resulting mesh
+    Info << "Writing resulting mesh ..." << endl;
     if (overwrite)
     {
         mesh.setInstance(oldInstance);
     }
     mesh.write();
 
-
     // Get list of cell splits.
     // (is for every cell in old mesh the cells they have been split into)
     const labelListList& oldToNew = multiRef.addedCells();
@@ -339,7 +380,7 @@ int main(int argc, char *argv[])
         }
     }
 
-    Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
+    Info<< "Writing refined cells (" << returnReduce(newCells.size(), sumOp<label>()) << ") to cellSet "
         << newCells.instance()/newCells.local()/newCells.name()
         << endl << endl;
 
@@ -390,7 +431,7 @@ int main(int argc, char *argv[])
         }
     }
 
-    Info<< "Writing map from new to old cell to "
+    Info<< "Writing cellMap (new to old): "
         << newToOld.objectPath() << nl << endl;
 
     newToOld.write();
-- 
1.9.1

diff --git a/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C b/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
index 493cd0c4a..7f6fc57a1 100644
--- a/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
+++ b/applications/utilities/mesh/manipulation/refineMesh/refineMesh.C
@@ -40,18 +40,11 @@ Description
 #include "argList.H"
 #include "polyMesh.H"
 #include "Time.H"
-#include "undoableMeshCutter.H"
-#include "hexCellLooper.H"
 #include "cellSet.H"
-#include "twoDPointCorrector.H"
-#include "directions.H"
-#include "OFstream.H"
 #include "multiDirRefinement.H"
 #include "labelIOList.H"
-#include "wedgePolyPatch.H"
-#include "plane.H"
-#include "SubField.H"
 #include "IOdictionary.H"
+#include "syncTools.H"
 
 using namespace Foam;
 
@@ -62,7 +55,7 @@ static const scalar edgeTol = 1e-3;
 
 
 // Print edge statistics on mesh.
-void printEdgeStats(const primitiveMesh& mesh)
+void printEdgeStats(const polyMesh& mesh)
 {
     label nX = 0;
     label nY = 0;
@@ -70,66 +63,90 @@ void printEdgeStats(const primitiveMesh& mesh)
 
     scalar minX = GREAT;
     scalar maxX = -GREAT;
-    vector x(1, 0, 0);
+    static const vector x(1, 0, 0);
 
     scalar minY = GREAT;
     scalar maxY = -GREAT;
-    vector y(0, 1, 0);
+    static const vector y(0, 1, 0);
 
     scalar minZ = GREAT;
     scalar maxZ = -GREAT;
-    vector z(0, 0, 1);
+    static const vector z(0, 0, 1);
 
     scalar minOther = GREAT;
     scalar maxOther = -GREAT;
 
+    PackedBoolList isMasterEdge(syncTools::getMasterEdges(mesh));
+
     const edgeList& edges = mesh.edges();
 
     forAll(edges, edgeI)
     {
-        const edge& e = edges[edgeI];
+        if (isMasterEdge[edgeI])
+        {
+            const edge& e = edges[edgeI];
 
-        vector eVec(e.vec(mesh.points()));
+            vector eVec(e.vec(mesh.points()));
 
-        scalar eMag = mag(eVec);
+            scalar eMag = mag(eVec);
 
-        eVec /= eMag;
+            eVec /= eMag;
 
-        if (mag(eVec & x) > 1-edgeTol)
-        {
-            minX = min(minX, eMag);
-            maxX = max(maxX, eMag);
-            nX++;
-        }
-        else if (mag(eVec & y) > 1-edgeTol)
-        {
-            minY = min(minY, eMag);
-            maxY = max(maxY, eMag);
-            nY++;
-        }
-        else if (mag(eVec & z) > 1-edgeTol)
-        {
-            minZ = min(minZ, eMag);
-            maxZ = max(maxZ, eMag);
-            nZ++;
-        }
-        else
-        {
-            minOther = min(minOther, eMag);
-            maxOther = max(maxOther, eMag);
+            if (mag(eVec & x) > 1-edgeTol)
+            {
+                minX = min(minX, eMag);
+                maxX = max(maxX, eMag);
+                nX++;
+            }
+            else if (mag(eVec & y) > 1-edgeTol)
+            {
+                minY = min(minY, eMag);
+                maxY = max(maxY, eMag);
+                nY++;
+            }
+            else if (mag(eVec & z) > 1-edgeTol)
+            {
+                minZ = min(minZ, eMag);
+                maxZ = max(maxZ, eMag);
+                nZ++;
+            }
+            else
+            {
+                minOther = min(minOther, eMag);
+                maxOther = max(maxOther, eMag);
+            }
         }
     }
 
-    Pout<< "Mesh edge statistics:" << endl
+    label nEdges = mesh.nEdges();
+    reduce(nEdges, sumOp<label>());
+    reduce(nX, sumOp<label>());
+    reduce(nY, sumOp<label>());
+    reduce(nZ, sumOp<label>());
+
+    reduce(minX, minOp<scalar>());
+    reduce(maxX, maxOp<scalar>());
+
+    reduce(minY, minOp<scalar>());
+    reduce(maxY, maxOp<scalar>());
+
+    reduce(minZ, minOp<scalar>());
+    reduce(maxZ, maxOp<scalar>());
+
+    reduce(minOther, minOp<scalar>());
+    reduce(maxOther, maxOp<scalar>());
+
+
+    Info<< "Mesh edge statistics:" << nl
         << "    x aligned :  number:" << nX << "\tminLen:" << minX
-        << "\tmaxLen:" << maxX << endl
+        << "\tmaxLen:" << maxX << nl
         << "    y aligned :  number:" << nY << "\tminLen:" << minY
-        << "\tmaxLen:" << maxY << endl
+        << "\tmaxLen:" << maxY << nl
         << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
-        << "\tmaxLen:" << maxZ << endl
-        << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
+        << "\tmaxLen:" << maxZ << nl
+        << "    other     :  number:" << nEdges - nX - nY - nZ
         << "\tminLen:" << minOther
-        << "\tmaxLen:" << maxOther << endl << endl;
+        << "\tmaxLen:" << maxOther << nl << endl;
 }
 
 
@@ -228,7 +245,8 @@ int main(int argc, char *argv[])
 
         cellSet cells(mesh, setName);
 
-        Pout<< "Read " << cells.size() << " cells from cellSet "
+        Info<< "Read " << returnReduce(cells.size(), sumOp<label>())
+            << " cells from cellSet "
             << cells.instance()/cells.local()/cells.name()
             << endl << endl;
 
@@ -239,12 +257,7 @@ int main(int argc, char *argv[])
         Info<< "Refining all cells" << nl << endl;
 
         // Select all cells
-        refCells.setSize(mesh.nCells());
-
-        forAll(mesh.cells(), celli)
-        {
-            refCells[celli] = celli;
-        }
+        refCells = identity(mesh.nCells());
 
         if (mesh.nGeometricD() == 3)
         {
@@ -340,7 +353,9 @@ int main(int argc, char *argv[])
         }
     }
 
-    Pout<< "Writing refined cells (" << newCells.size() << ") to cellSet "
+    Info<< "Writing refined cells ("
+        << returnReduce(newCells.size(), sumOp<label>())
+        << ") to cellSet "
         << newCells.instance()/newCells.local()/newCells.name()
         << endl << endl;
 
@@ -348,7 +363,6 @@ int main(int argc, char *argv[])
 
 
 
-
     //
     // Invert cell split to construct map from new to old
     //
diff --git a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
index d9aa92fdc..4ed3ff42d 100644
--- a/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
+++ b/src/dynamicMesh/meshCut/cellCuts/cellCuts.C
@@ -33,12 +33,25 @@ License
 #include "geomCellLooper.H"
 #include "OFstream.H"
 #include "plane.H"
+#include "syncTools.H"
+#include "dummyTransform.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
     defineTypeNameAndDebug(cellCuts, 0);
+
+    //- Template specialization for pTraits<edge> so we can use syncTools
+    //  functionality
+    template<>
+    class pTraits<edge>
+    {
+    public:
+
+        //- Component type
+        typedef edge cmptType;
+    };
 }
 
 
@@ -114,6 +127,148 @@ Foam::label Foam::cellCuts::firstUnique
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
+void Foam::cellCuts::syncProc()
+{
+    if (!Pstream::parRun())
+    {
+        return;
+    }
+
+    syncTools::syncPointList(mesh(), pointIsCut_, orEqOp<bool>(), false);
+    syncTools::syncEdgeList(mesh(), edgeIsCut_, orEqOp<bool>(), false);
+    syncTools::syncEdgeList(mesh(), edgeWeight_, maxEqOp<scalar>(), -GREAT);
+
+    {
+        const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
+
+        // Convert faceSplitCut into face-local data: vertex and edge w.r.t.
+        // vertex 0: (since this is same on both sides)
+        //
+        //      Sending-side vertex  Receiving-side vertex
+        //      0                   0
+        //      1                   3
+        //      2                   2
+        //      3                   1
+        //
+        //      Sending-side edge    Receiving side edge
+        //      0-1                  3-0
+        //      1-2                  2-3
+        //      2-3                  1-2
+        //      3-0                  0-1
+        //
+        // Encoding is as index:
+        //      0  : not set
+        //      >0 : vertex, vertex is index-1
+        //      <0 : edge, edge is -index-1
+
+        edgeList relCuts(nBnd, edge(0, 0));
+
+        const polyBoundaryMesh& pbm = mesh().boundaryMesh();
+
+        forAll(pbm, patchi)
+        {
+            const polyPatch& pp = pbm[patchi];
+
+            if (pp.coupled())
+            {
+                forAll(pp, i)
+                {
+                    label facei = pp.start()+i;
+                    label bFacei = facei-mesh().nInternalFaces();
+
+                    const Map<edge>::const_iterator iter =
+                        faceSplitCut_.find(facei);
+                    if (iter != faceSplitCut_.end())
+                    {
+                        const face& f = mesh().faces()[facei];
+                        const labelList& fEdges = mesh().faceEdges()[facei];
+                        const edge& cuts = iter();
+
+                        forAll(cuts, i)
+                        {
+                            if (isEdge(cuts[i]))
+                            {
+                                label edgei = getEdge(cuts[i]);
+                                label index = findIndex(fEdges, edgei);
+                                relCuts[bFacei][i] = -index-1;
+                            }
+                            else
+                            {
+                                label index = findIndex(f, getVertex(cuts[i]));
+                                relCuts[bFacei][i] = index+1;
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        // Exchange
+        syncTools::syncBoundaryFaceList
+        (
+            mesh(),
+            relCuts,
+            eqOp<edge>(),
+            dummyTransform()
+        );
+
+        // Convert relCuts back into mesh based data
+        forAll(pbm, patchi)
+        {
+            const polyPatch& pp = pbm[patchi];
+
+            if (pp.coupled())
+            {
+                forAll(pp, i)
+                {
+                    label facei = pp.start()+i;
+                    label bFacei = facei-mesh().nInternalFaces();
+
+                    const edge& relCut = relCuts[bFacei];
+                    if (relCut != edge(0, 0))
+                    {
+                        const face& f = mesh().faces()[facei];
+                        const labelList& fEdges = mesh().faceEdges()[facei];
+
+                        edge absoluteCut(0, 0);
+                        forAll(relCut, i)
+                        {
+                            if (relCut[i] < 0)
+                            {
+                                label oppFp = -relCut[i]-1;
+                                label fp = f.size()-1-oppFp;
+                                absoluteCut[i] = edgeToEVert(fEdges[fp]);
+                            }
+                            else
+                            {
+                                label oppFp = relCut[i]-1;
+                                label fp = f.size()-1-oppFp;
+                                absoluteCut[i] = vertToEVert(f[fp]);
+                            }
+                        }
+
+                        if
+                        (
+                           !faceSplitCut_.insert(facei, absoluteCut)
+                         && faceSplitCut_[facei] != absoluteCut
+                        )
+                        {
+                            FatalErrorInFunction
+                                << "Cut " << faceSplitCut_[facei]
+                                << " on face " << mesh().faceCentres()[facei]
+                                << " of coupled patch " << pp.name()
+                                << " is not consistent with coupled cut "
+                                << absoluteCut
+                                << exit(FatalError);
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::cellCuts::writeUncutOBJ
 (
     const fileName& dir,
@@ -337,7 +492,7 @@ Foam::label Foam::cellCuts::vertexVertexToFace
 
 void Foam::cellCuts::calcFaceCuts() const
 {
-    if (faceCutsPtr_)
+    if (faceCutsPtr_.valid())
     {
         FatalErrorInFunction
             << "faceCuts already calculated" << abort(FatalError);
@@ -345,9 +500,8 @@ void Foam::cellCuts::calcFaceCuts() const
 
     const faceList& faces = mesh().faces();
 
-
-    faceCutsPtr_ = new labelListList(mesh().nFaces());
-    labelListList& faceCuts = *faceCutsPtr_;
+    faceCutsPtr_.reset(new labelListList(mesh().nFaces()));
+    labelListList& faceCuts = faceCutsPtr_();
 
     for (label facei = 0; facei < mesh().nFaces(); facei++)
     {
@@ -2624,6 +2778,16 @@ void Foam::cellCuts::check() const
 
 
     // Check that cut faces have a neighbour that is cut.
+    boolList nbrCellIsCut;
+    {
+        boolList cellIsCut(mesh().nCells(), false);
+        forAll(cellLoops_, celli)
+        {
+            cellIsCut[celli] = cellLoops_[celli].size();
+        }
+        syncTools::swapBoundaryCellList(mesh(), cellIsCut, nbrCellIsCut);
+    }
+
     forAllConstIter(Map<edge>, faceSplitCut_, iter)
     {
         label facei = iter.key();
@@ -2645,9 +2809,11 @@ void Foam::cellCuts::check() const
         }
         else
         {
+            label bFacei = facei - mesh().nInternalFaces();
+
             label own = mesh().faceOwner()[facei];
 
-            if (cellLoops_[own].empty())
+            if (cellLoops_[own].empty() && !nbrCellIsCut[bFacei])
             {
                 FatalErrorInFunction
                     << "Boundary face:" << facei << " cut by " << iter()
@@ -2675,7 +2841,6 @@ Foam::cellCuts::cellCuts
     pointIsCut_(expand(mesh.nPoints(), meshVerts)),
     edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
     edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
-    faceCutsPtr_(nullptr),
     faceSplitCut_(cutCells.size()),
     cellLoops_(mesh.nCells()),
     nLoops_(-1),
@@ -2718,7 +2883,6 @@ Foam::cellCuts::cellCuts
     pointIsCut_(expand(mesh.nPoints(), meshVerts)),
     edgeIsCut_(expand(mesh.nEdges(), meshEdges)),
     edgeWeight_(expand(mesh.nEdges(), meshEdges, meshEdgeWeights)),
-    faceCutsPtr_(nullptr),
     faceSplitCut_(mesh.nFaces()/10 + 1),
     cellLoops_(mesh.nCells()),
     nLoops_(-1),
@@ -2734,6 +2898,9 @@ Foam::cellCuts::cellCuts
 
     calcLoopsAndAddressing(identity(mesh.nCells()));
 
+    // Adds cuts on other side of coupled boundaries
+    syncProc();
+
     // Calculate planes and flip cellLoops if necessary
     orientPlanesAndLoops();
 
@@ -2763,7 +2930,6 @@ Foam::cellCuts::cellCuts
     pointIsCut_(mesh.nPoints(), false),
     edgeIsCut_(mesh.nEdges(), false),
     edgeWeight_(mesh.nEdges(), -GREAT),
-    faceCutsPtr_(nullptr),
     faceSplitCut_(cellLabels.size()),
     cellLoops_(mesh.nCells()),
     nLoops_(-1),
@@ -2778,6 +2944,9 @@ Foam::cellCuts::cellCuts
     // Makes sure cuts are consistent
     setFromCellLoops(cellLabels, cellLoops, cellEdgeWeights);
 
+    // Adds cuts on other side of coupled boundaries
+    syncProc();
+
     // Calculate planes and flip cellLoops if necessary
     orientPlanesAndLoops();
 
@@ -2806,7 +2975,6 @@ Foam::cellCuts::cellCuts
     pointIsCut_(mesh.nPoints(), false),
     edgeIsCut_(mesh.nEdges(), false),
     edgeWeight_(mesh.nEdges(), -GREAT),
-    faceCutsPtr_(nullptr),
     faceSplitCut_(refCells.size()),
     cellLoops_(mesh.nCells()),
     nLoops_(-1),
@@ -2821,6 +2989,9 @@ Foam::cellCuts::cellCuts
     // Makes sure cuts are consistent
     setFromCellCutter(cellCutter, refCells);
 
+    // Adds cuts on other side of coupled boundaries
+    syncProc();
+
     // Calculate planes and flip cellLoops if necessary
     orientPlanesAndLoops();
 
@@ -2850,7 +3021,6 @@ Foam::cellCuts::cellCuts
     pointIsCut_(mesh.nPoints(), false),
     edgeIsCut_(mesh.nEdges(), false),
     edgeWeight_(mesh.nEdges(), -GREAT),
-    faceCutsPtr_(nullptr),
     faceSplitCut_(cellLabels.size()),
     cellLoops_(mesh.nCells()),
     nLoops_(-1),
@@ -2866,6 +3036,9 @@ Foam::cellCuts::cellCuts
     // Makes sure cuts are consistent
     setFromCellCutter(cellCutter, cellLabels, cutPlanes);
 
+    // Adds cuts on other side of coupled boundaries
+    syncProc();
+
     // Calculate planes and flip cellLoops if necessary
     orientPlanesAndLoops();
 
@@ -2900,7 +3073,6 @@ Foam::cellCuts::cellCuts
     pointIsCut_(pointIsCut),
     edgeIsCut_(edgeIsCut),
     edgeWeight_(edgeWeight),
-    faceCutsPtr_(nullptr),
     faceSplitCut_(faceSplitCut),
     cellLoops_(cellLoops),
     nLoops_(nLoops),
@@ -2924,7 +3096,7 @@ Foam::cellCuts::~cellCuts()
 
 void Foam::cellCuts::clearOut()
 {
-    deleteDemandDrivenData(faceCutsPtr_);
+    faceCutsPtr_.clear();
 }
 
 
diff --git a/src/dynamicMesh/meshCut/cellCuts/cellCuts.H b/src/dynamicMesh/meshCut/cellCuts/cellCuts.H
index 7e80ce18f..fed3b9308 100644
--- a/src/dynamicMesh/meshCut/cellCuts/cellCuts.H
+++ b/src/dynamicMesh/meshCut/cellCuts/cellCuts.H
@@ -128,7 +128,7 @@ class cellCuts
 
             //- Cuts per existing face (includes those along edge of face)
             //  Cuts in no particular order.
-            mutable labelListList* faceCutsPtr_;
+            mutable autoPtr<labelListList> faceCutsPtr_;
 
             //- Per face : cut across edge (so not along existing edge)
             //  (can only be one per face)
@@ -418,6 +418,10 @@ class cellCuts
                 const List<scalarField>& cellLoopWeights
             );
 
+            //- Update basic cut information to be consistent across
+            //  coupled boundaries.
+            void syncProc();
+
             //- Cut cells and update basic cut information from cellLoops.
             //  Checks each loop for consistency with existing cut pattern.
             void setFromCellCutter
@@ -557,11 +561,11 @@ public:
             //  Cuts in no particular order
             const labelListList& faceCuts() const
             {
-                if (!faceCutsPtr_)
+                if (!faceCutsPtr_.valid())
                 {
                     calcFaceCuts();
                 }
-                return *faceCutsPtr_;
+                return faceCutsPtr_();
             }
 
             //- Gives for split face the two cuts that split the face into two.
diff --git a/src/dynamicMesh/meshCut/directions/directions.C b/src/dynamicMesh/meshCut/directions/directions.C
index 2097a2c57..4dd64f62e 100644
--- a/src/dynamicMesh/meshCut/directions/directions.C
+++ b/src/dynamicMesh/meshCut/directions/directions.C
@@ -257,7 +257,11 @@ Foam::vectorField Foam::directions::propagateDirection
         }
     }
 
-    Pout<< "Calculated local coords for " << defaultDir
+    reduce(nGeom, sumOp<label>());
+    reduce(nTopo, sumOp<label>());
+    reduce(nUnset, sumOp<label>());
+
+    Info<< "Calculated local coords for " << defaultDir
         << endl
         << "    Geometric cut cells   : " << nGeom << endl
         << "    Topological cut cells : " << nTopo << endl
@@ -322,7 +326,7 @@ Foam::directions::directions
         vector normal = tan1 ^ tan2;
         normal /= mag(normal);
 
-        Pout<< "Global Coordinate system:" << endl
+        Info<< "Global Coordinate system:" << endl
             << "     normal : " << normal << endl
             << "     tan1   : " << tan1 << endl
             << "     tan2   : " << tan2
diff --git a/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C b/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C
index 903e1c60a..025157c3a 100644
--- a/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C
+++ b/src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C
@@ -33,6 +33,7 @@ License
 #include "polyAddPoint.H"
 #include "polyAddFace.H"
 #include "polyAddCell.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -536,7 +537,7 @@ void Foam::meshCutter::setRefinement
     addedPoints_.clear();
     addedPoints_.resize(cuts.nLoops());
 
-    if (cuts.nLoops() == 0)
+    if (returnReduce(cuts.nLoops(), sumOp<label>()) == 0)
     {
         return;
     }
@@ -544,6 +545,40 @@ void Foam::meshCutter::setRefinement
     const labelListList& anchorPts = cuts.cellAnchorPoints();
     const labelListList& cellLoops = cuts.cellLoops();
 
+    if (debug)
+    {
+        // Check that any edge is cut only if any cell using it is cut
+        boolList edgeOnCutCell(mesh().nEdges(), false);
+        forAll(cuts.cellLoops(), celli)
+        {
+            if (cuts.cellLoops()[celli].size())
+            {
+                const labelList& cEdges = mesh().cellEdges(celli);
+                forAll(cEdges, i)
+                {
+                    edgeOnCutCell[cEdges[i]] = true;
+                }
+            }
+        }
+        syncTools::syncEdgeList(mesh(), edgeOnCutCell, orEqOp<bool>(), false);
+
+        forAll(cuts.edgeIsCut(), edgeI)
+        {
+            if (cuts.edgeIsCut()[edgeI] && !edgeOnCutCell[edgeI])
+            {
+                const edge& e = mesh().edges()[edgeI];
+
+                WarningInFunction
+                    << "Problem: cut edge but none of the cells using"
+                    << " it is cut\n"
+                    << "edge:" << edgeI << " verts:" << e
+                    << " at:" << e.line(mesh().points())
+                    << endl;    //abort(FatalError);
+            }
+        }
+    }
+
+
     //
     // Add new points along cut edges.
     //
@@ -554,15 +589,6 @@ void Foam::meshCutter::setRefinement
         {
             const edge& e = mesh().edges()[edgeI];
 
-            // Check if there is any cell using this edge.
-            if (debug && findCutCell(cuts, mesh().edgeCells()[edgeI]) == -1)
-            {
-                FatalErrorInFunction
-                    << "Problem: cut edge but none of the cells using it is\n"
-                    << "edge:" << edgeI << " verts:" << e
-                    << abort(FatalError);
-            }
-
             // One of the edge end points should be master point of nbCelli.
             label masterPointi = e.start();
 
diff --git a/src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C b/src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C
index 27ee747ae..c5d1b5526 100644
--- a/src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C
+++ b/src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C
@@ -270,7 +270,7 @@ Foam::Map<Foam::label> Foam::refinementIterator::setRefinement
     while (!stop);
 
 
-    if (nRefCells == oldRefCells)
+    if (returnReduce((nRefCells == oldRefCells), andOp<bool>()))
     {
         WarningInFunction
             << "stopped refining."