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