2017-12-15 17:58 GMT

View Issue Details Jump to Notes ]
IDProjectCategoryView StatusLast Update
0002621OpenFOAM[All Projects] Bugpublic2017-08-15 13:10
Reportersose 
Assigned Tohenry 
PrioritynormalSeveritycrashReproducibilityalways
StatusresolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version14.04
Product Version 
Target VersionFixed in Versiondev 
Summary0002621: refineMesh crashes when run in parallel
DescriptionrefineMesh does not do parallel communication correctly causing non-matching processor patches after a cut. To demonstrate this behavior please see the attached test case. Allrun will run "blockMesh", "setSet" and "refineMesh" successfully in serial, but fail in parallel.
Steps To Reproduceexecute Allrun from the attached test case
TagsNo tags attached.
Attached Files
  • tgz file icon bug-refineMesh-in-parallel.tgz (2,695 bytes) 2017-07-19 09:05
  • patch file icon BUG2621-refineMesh-in-parallel.patch (14,471 bytes) 2017-07-19 10:48 -
    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
    
    
    patch file icon BUG2621-refineMesh-in-parallel.patch (14,471 bytes) 2017-07-19 10:48 +
  • patch file icon refineMesh.patch (22,618 bytes) 2017-07-21 17:09 -
    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."
    
    patch file icon refineMesh.patch (22,618 bytes) 2017-07-21 17:09 +

-Relationships
+Relationships

-Notes

~0008412

henry (manager)

Can you provide a patch to fix the problem?

~0008413

sose (reporter)

The problem seems to be in the cellCuts class. I add a few lines of code to sync. its private data. Please see the attachment. The problem seems to disappear ... at least for this specific case where this constructor is used

        cellCuts
        (
            const polyMesh& mesh,
            const cellLooper& cellCutter,
            const List<refineCell>& refCells
        );

~0008431

MattijsJ (reporter)

Impressive work. I could not apply your patch but I think I've done something similar. Could you try this one?

~0008451

sose (reporter)

Obviously you have put more thought into fixing this bug than I do. When you check whether or not a patch is a processor patch I think it is better to use directly isA<processorPolyPatch> rather than "coupled()". If I understand correctly "coupled()" might not be "processorPolyPatch". Anyway, I've tested your patch and it works fine. Thank you for sharing.

~0008590

henry (manager)

Resolved by commits

c6f365f6ccc401e3d4a45ff6af332a7673bb7a3c
78e79c4ad309de969df66571488a123d3a5064fb
+Notes

-Issue History
Date Modified Username Field Change
2017-07-19 09:05 sose New Issue
2017-07-19 09:05 sose File Added: bug-refineMesh-in-parallel.tgz
2017-07-19 09:14 henry Note Added: 0008412
2017-07-19 10:48 sose File Added: BUG2621-refineMesh-in-parallel.patch
2017-07-19 10:48 sose Note Added: 0008413
2017-07-21 17:09 MattijsJ File Added: refineMesh.patch
2017-07-21 17:09 MattijsJ Note Added: 0008431
2017-07-24 10:39 sose Note Added: 0008451
2017-08-15 13:10 henry Assigned To => henry
2017-08-15 13:10 henry Status new => resolved
2017-08-15 13:10 henry Resolution open => fixed
2017-08-15 13:10 henry Fixed in Version => dev
2017-08-15 13:10 henry Note Added: 0008590
+Issue History