refineMesh: Correct parallel operation

Patch contributed by Mattijs Janssens
Resolves bug-report https://bugs.openfoam.org/view.php?id=2621

applications/utilities/mesh/manipulation/refineMesh/refineMesh.C

@@ -43,18 +43,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;
 
@@ -65,7 +58,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;
@@ -73,22 +66,26 @@ 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)
+    {
+        if (isMasterEdge[edgeI])
         {
             const edge& e = edges[edgeI];
 
@@ -122,17 +119,37 @@ void printEdgeStats(const primitiveMesh& mesh)
                 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
+        << "\tmaxLen:" << maxZ << nl
-        << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
+        << "    other     :  number:" << nEdges - nX - nY - nZ
         << "\tminLen:" << minOther
-        << "\tmaxLen:" << maxOther << endl << endl;
+        << "\tmaxLen:" << maxOther << nl << endl;
 }
 
 
@@ -231,7 +248,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;
 
@@ -242,12 +260,7 @@ int main(int argc, char *argv[])
         Info<< "Refining all cells" << nl << endl;
 
         // Select all cells
-        refCells.setSize(mesh.nCells());
+        refCells = identity(mesh.nCells());
-
-        forAll(mesh.cells(), celli)
-        {
-            refCells[celli] = celli;
-        }
 
         if (mesh.nGeometricD() == 3)
         {
@@ -343,7 +356,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;
 
@@ -351,7 +366,6 @@ int main(int argc, char *argv[])
 
 
 
-
     //
     // Invert cell split to construct map from new to old
     //

src/dynamicMesh/meshCut/cellCuts/cellCuts.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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_.reset(new labelListList(mesh().nFaces()));
-    faceCutsPtr_ = new labelListList(mesh().nFaces());
+    labelListList& faceCuts = faceCutsPtr_();
-    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();
 }
 
 

src/dynamicMesh/meshCut/cellCuts/cellCuts.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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.

src/dynamicMesh/meshCut/directions/directions.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -250,7 +250,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
@@ -315,7 +319,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

src/dynamicMesh/meshCut/meshModifiers/meshCutter/meshCutter.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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();
 

src/dynamicMesh/meshCut/meshModifiers/refinementIterator/refinementIterator.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -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."