ENH: foamyHexMesh: cell/face zone in parallel

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -600,6 +600,19 @@ private:
             PackedBoolList& boundaryFacesToRemove
         );
 
+        void calcNeighbourCellCentres
+        (
+            const polyMesh& mesh,
+            const pointField& cellCentres,
+            pointField& neiCc
+        ) const;
+
+        void selectSeparatedCoupledFaces
+        (
+            const polyMesh& mesh,
+            boolList& selected
+        ) const;
+
         //- From meshRefinementBaffles.C. Use insidePoint for a surface to
         //  determine the cell zone.
         void findCellZoneInsideWalk

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -2744,12 +2744,45 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                 }
                 else
                 {
-                    // internal face
+//                    if
-                    faces[dualFaceI] = newDualFace;
+//                    (
-                    owner[dualFaceI] = own;
+//                        ptPairs_.isPointPair(vA, vB)
-                    neighbour[dualFaceI] = nei;
+//                     || ftPtConformer_.featurePointPairs().isPointPair(vA, vB)
+//                    )
+//                    {
+                        patchIndex = geometryToConformTo_.findPatch(ptA, ptB);
+//                    }
 
-                    dualFaceI++;
+                    if (patchIndex != -1)
+                    {
+//                        patchFaces[patchIndex].append(newDualFace);
+//                        patchOwners[patchIndex].append(own);
+//                        indirectPatchFace[patchIndex].append(false);
+//
+//                        if
+//                        (
+//                            labelPair(vB->index(), vB->procIndex())
+//                          < labelPair(vA->index(), vA->procIndex())
+//                        )
+//                        {
+//                            patchPPSlaves[patchIndex].append(vB->index());
+//                        }
+//                        else
+//                        {
+//                            patchPPSlaves[patchIndex].append(vA->index());
+//                        }
+
+//                        baffleFaces[dualFaceI] = patchIndex;
+                    }
+//                    else
+                    {
+                        // internal face
+                        faces[dualFaceI] = newDualFace;
+                        owner[dualFaceI] = own;
+                        neighbour[dualFaceI] = nei;
+
+                        dualFaceI++;
+                    }
                 }
             }
         }

applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C

@@ -111,7 +111,7 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
 
     {
         pointField points;
-        labelList boundaryPts(number_of_finite_cells(), -1);
+        labelList boundaryPts;
         faceList faces;
         labelList owner;
         labelList neighbour;
@@ -407,6 +407,78 @@ void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
 }
 
 
+void Foam::conformalVoronoiMesh::calcNeighbourCellCentres
+(
+    const polyMesh& mesh,
+    const pointField& cellCentres,
+    pointField& neiCc
+) const
+{
+    label nBoundaryFaces = mesh.nFaces() - mesh.nInternalFaces();
+
+    if (neiCc.size() != nBoundaryFaces)
+    {
+        FatalErrorIn("conformalVoronoiMesh::calcNeighbourCellCentres(..)")
+            << "nBoundaries:" << nBoundaryFaces
+            << " neiCc:" << neiCc.size()
+            << abort(FatalError);
+    }
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        const labelUList& faceCells = pp.faceCells();
+
+        label bFaceI = pp.start() - mesh.nInternalFaces();
+
+        if (pp.coupled())
+        {
+            forAll(faceCells, i)
+            {
+                neiCc[bFaceI] = cellCentres[faceCells[i]];
+                bFaceI++;
+            }
+        }
+    }
+
+    // Swap coupled boundaries. Apply separation to cc since is coordinate.
+    syncTools::swapBoundaryFacePositions(mesh, neiCc);
+}
+
+
+void Foam::conformalVoronoiMesh::selectSeparatedCoupledFaces
+(
+    const polyMesh& mesh,
+    boolList& selected
+) const
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        // Check all coupled. Avoid using .coupled() so we also pick up AMI.
+        if (isA<coupledPolyPatch>(patches[patchI]))
+        {
+            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+            (
+                patches[patchI]
+            );
+
+            if (cpp.separated() || !cpp.parallel())
+            {
+                forAll(cpp, i)
+                {
+                    selected[cpp.start()+i] = true;
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::conformalVoronoiMesh::findCellZoneInsideWalk
 (
     const polyMesh& mesh,
@@ -417,7 +489,7 @@ void Foam::conformalVoronoiMesh::findCellZoneInsideWalk
 {
     // Analyse regions. Reuse regionsplit
     boolList blockedFace(mesh.nFaces());
-    //selectSeparatedCoupledFaces(blockedFace);
+    selectSeparatedCoupledFaces(mesh, blockedFace);
 
     forAll(faceToSurface, faceI)
     {
@@ -630,42 +702,90 @@ void Foam::conformalVoronoiMesh::calcFaceZones
     const labelList& faceOwner = mesh.faceOwner();
     const labelList& faceNeighbour = mesh.faceNeighbour();
 
+    labelList neiFaceOwner(mesh.nFaces() - mesh.nInternalFaces(), -1);
+
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        const labelUList& faceCells = pp.faceCells();
+
+        label bFaceI = pp.start() - mesh.nInternalFaces();
+
+        if (pp.coupled())
+        {
+            forAll(faceCells, i)
+            {
+                neiFaceOwner[bFaceI] = cellToSurface[faceCells[i]];
+                bFaceI++;
+            }
+        }
+    }
+
+    syncTools::swapBoundaryFaceList(mesh, neiFaceOwner);
+
     forAll(faces, faceI)
     {
         const label ownerSurfaceI = cellToSurface[faceOwner[faceI]];
 
+        if (faceToSurface[faceI] >= 0)
+        {
+            continue;
+        }
+
         if (mesh.isInternalFace(faceI))
         {
             const label neiSurfaceI = cellToSurface[faceNeighbour[faceI]];
 
-            flipMap[faceI] =
-                (
-                    ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
-                  ? false
-                  : true
-                );
-
             if
             (
                 (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
              && ownerSurfaceI != neiSurfaceI
             )
             {
-                if (ownerSurfaceI > neiSurfaceI)
+                flipMap[faceI] =
-                {
+                    (
-                    faceToSurface[faceI] = ownerSurfaceI;
+                        ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
-                }
+                      ? false
-                else
+                      : true
-                {
+                    );
-                    faceToSurface[faceI] = neiSurfaceI;
+
-                }
+                faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
             }
         }
         else
         {
-            if (ownerSurfaceI >= 0)
+            label patchID = mesh.boundaryMesh().whichPatch(faceI);
+
+            if (mesh.boundaryMesh()[patchID].coupled())
             {
-                faceToSurface[faceI] = ownerSurfaceI;
+                const label neiSurfaceI =
+                    neiFaceOwner[faceI - mesh.nInternalFaces()];
+
+                if
+                (
+                    (ownerSurfaceI >= 0 || neiSurfaceI >= 0)
+                 && ownerSurfaceI != neiSurfaceI
+                )
+                {
+                    flipMap[faceI] =
+                        (
+                            ownerSurfaceI == max(ownerSurfaceI, neiSurfaceI)
+                          ? false
+                          : true
+                        );
+
+                    faceToSurface[faceI] = max(ownerSurfaceI, neiSurfaceI);
+                }
+            }
+            else
+            {
+                if (ownerSurfaceI >= 0)
+                {
+                    faceToSurface[faceI] = ownerSurfaceI;
+                }
             }
         }
     }
@@ -674,104 +794,150 @@ void Foam::conformalVoronoiMesh::calcFaceZones
     const PtrList<surfaceZonesInfo>& surfZones =
         geometryToConformTo().surfZones();
 
-    labelList insidePointNamedSurfaces
+    labelList unclosedSurfaces
     (
-        surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
+        surfaceZonesInfo::getUnclosedNamedSurfaces
+        (
+            surfZones,
+            geometryToConformTo().geometry(),
+            geometryToConformTo().surfaces()
+        )
     );
 
+    pointField neiCc(mesh.nFaces() - mesh.nInternalFaces());
+    calcNeighbourCellCentres
+    (
+        mesh,
+        cellCentres,
+        neiCc
+    );
+
+    OBJstream intersections(time().path()/"ints.obj");
+    OBJstream intersectionFaces(time().path()/"intFaces.obj");
+
     // Use intersection of cellCentre connections
     forAll(faces, faceI)
     {
-        if
+        if (faceToSurface[faceI] >= 0)
-        (
+        {
-            mesh.isInternalFace(faceI)
+            continue;
-         && faceToSurface[faceI] < 0
+        }
-        )
+
+        label patchID = mesh.boundaryMesh().whichPatch(faceI);
+
+        const label own = faceOwner[faceI];
+
+        List<pointIndexHit> surfHit;
+        labelList hitSurface;
+
+        if (mesh.isInternalFace(faceI))
         {
-            const label own = faceOwner[faceI];
             const label nei = faceNeighbour[faceI];
 
-            pointIndexHit surfHit;
+            geometryToConformTo().findSurfaceAllIntersections
-            label hitSurface;
-
-            geometryToConformTo().findSurfaceAnyIntersection
             (
                 cellCentres[own],
                 cellCentres[nei],
                 surfHit,
                 hitSurface
             );
+        }
+        else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
+        {
+            geometryToConformTo().findSurfaceAllIntersections
+            (
+                cellCentres[own],
+                neiCc[faceI - mesh.nInternalFaces()],
+                surfHit,
+                hitSurface
+            );
 
-            if (surfHit.hit())
+            if (surfHit.size() == 1 && surfHit[0].hit())
             {
-                if (findIndex(insidePointNamedSurfaces, hitSurface) != -1)
+                intersections.write
-                {
+                (
-                    faceToSurface[faceI] = hitSurface;
+                    linePointRef
-
-                    vectorField norm;
-                    geometryToConformTo().getNormal
                     (
-                        hitSurface,
+                        cellCentres[own],
-                        List<pointIndexHit>(1, surfHit),
+                        neiCc[faceI - mesh.nInternalFaces()]
-                        norm
+                    )
-                    );
+                );
-
-                    const vector fN = faces[faceI].normal(mesh.points());
-
-                    if ((norm[0] & fN/(mag(fN) + SMALL)) < 0)
-                    {
-                        flipMap[faceI] = true;
-                    }
-                    else
-                    {
-                        flipMap[faceI] = false;
-                    }
-                }
             }
         }
-    }
 
-
+        // If there are multiple intersections then do not add to
-    labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
+        // a faceZone
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+        if (surfHit.size() == 1 && surfHit[0].hit())
-
-    forAll(patches, patchI)
-    {
-        const polyPatch& pp = patches[patchI];
-
-        if (pp.coupled())
         {
-            forAll(pp, i)
+            if (findIndex(unclosedSurfaces, hitSurface[0]) != -1)
             {
-                label faceI = pp.start()+i;
+                vectorField norm;
-                label ownSurface = cellToSurface[faceOwner[faceI]];
+                geometryToConformTo().getNormal
-                neiCellSurface[faceI - mesh.nInternalFaces()] = ownSurface;
+                (
-            }
+                    hitSurface[0],
-        }
+                    List<pointIndexHit>(1, surfHit[0]),
-    }
+                    norm
-    syncTools::swapBoundaryFaceList(mesh, neiCellSurface);
+                );
 
-    forAll(patches, patchI)
+                vector fN = faces[faceI].normal(mesh.points());
-    {
+                fN /= mag(fN) + SMALL;
-        const polyPatch& pp = patches[patchI];
 
-        if (pp.coupled())
+                if ((norm[0] & fN) < 0)
-        {
-            forAll(pp, i)
-            {
-                label faceI = pp.start()+i;
-                label ownSurface = cellToSurface[faceOwner[faceI]];
-                label neiSurface = neiCellSurface[faceI-mesh.nInternalFaces()];
-
-                if (faceToSurface[faceI] == -1 && (ownSurface != neiSurface))
                 {
-                    // Give face the max cell zone
+                    flipMap[faceI] = true;
-                    faceToSurface[faceI] =  max(ownSurface, neiSurface);
                 }
+                else
+                {
+                    flipMap[faceI] = false;
+                }
+
+                faceToSurface[faceI] = hitSurface[0];
+
+                intersectionFaces.write(faces[faceI], mesh.points());
             }
         }
     }
 
+
+//    labelList neiCellSurface(mesh.nFaces()-mesh.nInternalFaces());
+//
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (pp.coupled())
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//                label ownSurface = cellToSurface[faceOwner[faceI]];
+//                neiCellSurface[faceI - mesh.nInternalFaces()] = ownSurface;
+//            }
+//        }
+//    }
+//    syncTools::swapBoundaryFaceList(mesh, neiCellSurface);
+//
+//    forAll(patches, patchI)
+//    {
+//        const polyPatch& pp = patches[patchI];
+//
+//        if (pp.coupled())
+//        {
+//            forAll(pp, i)
+//            {
+//                label faceI = pp.start()+i;
+//                label ownSurface = cellToSurface[faceOwner[faceI]];
+//                label neiSurface = neiCellSurface[faceI-mesh.nInternalFaces()];
+//
+//                if (faceToSurface[faceI] == -1 && (ownSurface != neiSurface))
+//                {
+//                    // Give face the max cell zone
+//                    faceToSurface[faceI] =  max(ownSurface, neiSurface);
+//                }
+//            }
+//        }
+//    }
+//
     // Sync
     syncTools::syncFaceList(mesh, faceToSurface, maxEqOp<label>());
 }
@@ -1049,7 +1215,7 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
 
     pBufs.finishedSends();
 
-    Info<< incrIndent << indent << "    Face ordering initialised..." << endl;
+    Info<< incrIndent << indent << "Face ordering initialised..." << endl;
 
     // Receive and calculate ordering
     bool anyChanged = false;
@@ -1108,7 +1274,7 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
         }
     }
 
-    Info<< incrIndent << indent << "    Faces matched." << endl;
+    Info<< incrIndent << indent << "Faces matched." << endl;
 
     reduce(anyChanged, orOp<bool>());
 
@@ -1145,6 +1311,7 @@ void Foam::conformalVoronoiMesh::reorderProcessorPatches
             << " faces have been reordered" << nl
             << indent << returnReduce(nRotated, sumOp<label>())
             << " faces have been rotated"
+            << decrIndent << decrIndent
             << decrIndent << decrIndent << endl;
     }
 }
@@ -1192,7 +1359,8 @@ void Foam::conformalVoronoiMesh::writeMesh
         );
     }
 
-    Info<< "    Constructing mesh" << endl;
+    Info<< incrIndent;
+    Info<< indent << "Constructing mesh" << endl;
 
     timeCheck("Before fvMesh construction");
 
@@ -1212,7 +1380,7 @@ void Foam::conformalVoronoiMesh::writeMesh
         xferMove(neighbour)
     );
 
-    Info<< "    Adding patches to mesh" << endl;
+    Info<< indent << "Adding patches to mesh" << endl;
 
     List<polyPatch*> patches(patchNames.size());
 
@@ -1258,11 +1426,7 @@ void Foam::conformalVoronoiMesh::writeMesh
             // Check that the patch is not empty on every processor
             reduce(totalPatchSize, sumOp<label>());
 
-            if
+            if (totalPatchSize > 0)
-            (
-                totalPatchSize > 0
-//             && !geometryToConformTo().surfZones().set(p)
-            )
             {
                 patches[nValidPatches] = polyPatch::New
                 (
@@ -1377,16 +1541,40 @@ void Foam::conformalVoronoiMesh::writeMesh
         const labelList& faceOwner = mesh.faceOwner();
         const labelList& faceNeighbour = mesh.faceNeighbour();
 
+//        // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
+//        labelList neiCellZone(mesh.nFaces() - mesh.nInternalFaces(), -1);
+//
+//        const polyBoundaryMesh& patches = mesh.boundaryMesh();
+//
+//        forAll(patches, patchI)
+//        {
+//            const polyPatch& pp = patches[patchI];
+//
+//            if (pp.coupled())
+//            {
+//                forAll(pp, i)
+//                {
+//                    label faceI = pp.start()+i;
+//                    neiCellZone[faceI - mesh.nInternalFaces()] =
+//                        surfaceToCellZone[cellToSurface[faceOwner[faceI]]];
+//                }
+//            }
+//        }
+//
+//        syncTools::swapBoundaryFaceList(mesh, neiCellZone);
+
         forAll(faceToSurface, faceI)
         {
-            if (!mesh.isInternalFace(faceI))
+            label surfaceI = faceToSurface[faceI];
+
+            if (surfaceI < 0)
             {
                 continue;
             }
 
-            label surfaceI = faceToSurface[faceI];
+            label patchID = mesh.boundaryMesh().whichPatch(faceI);
 
-            if (surfaceI >= 0)
+            if (mesh.isInternalFace(faceI))
             {
                 label own = faceOwner[faceI];
                 label nei = faceNeighbour[faceI];
@@ -1407,6 +1595,26 @@ void Foam::conformalVoronoiMesh::writeMesh
                     )
                 );
             }
+            else if (patchID != -1 && mesh.boundaryMesh()[patchID].coupled())
+            {
+                label own = faceOwner[faceI];
+
+                meshMod.setAction
+                (
+                    polyModifyFace
+                    (
+                        mesh.faces()[faceI],           // modified face
+                        faceI,                          // label of face
+                        own,               // owner
+                        -1,                             // neighbour
+                        false,                          // face flip
+                        patchID,                         // patch for face
+                        false,                          // remove from zone
+                        surfaceToFaceZone[surfaceI],       // zone for face
+                        flipMap[faceI]                            // face flip in zone
+                    )
+                );
+            }
         }
 
         // Change the mesh (no inflation, parallel sync)

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalMeshData.C

@@ -847,7 +847,9 @@ Foam::label Foam::globalMeshData::findTransform
     {
         FatalErrorIn("globalMeshData::findTransform(..)")
             << "Problem. Cannot find " << remotePoint
-            << " or " << localPoint << " in " << info
+            << " or " << localPoint  << " "
+            << coupledPatch().localPoints()[localPoint]
+            << " in " << info
             << endl
             << "remoteTransformI:" << remoteTransformI << endl
             << "localTransformI:" << localTransformI

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C

@@ -182,7 +182,6 @@ Foam::surfaceZonesInfo::surfaceZonesInfo(const surfaceZonesInfo& surfZone)
 {}
 
 
-// Get indices of unnamed surfaces (surfaces without faceZoneName)
 Foam::labelList Foam::surfaceZonesInfo::getUnnamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList
@@ -204,7 +203,6 @@ Foam::labelList Foam::surfaceZonesInfo::getUnnamedSurfaces
 }
 
 
-// Get indices of named surfaces (surfaces with faceZoneName)
 Foam::labelList Foam::surfaceZonesInfo::getNamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList
@@ -230,7 +228,6 @@ Foam::labelList Foam::surfaceZonesInfo::getNamedSurfaces
 }
 
 
-// Get indices of closed named surfaces
 Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList,
@@ -263,7 +260,33 @@ Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
 }
 
 
-// Get indices of closed named surfaces
+Foam::labelList Foam::surfaceZonesInfo::getUnclosedNamedSurfaces
+(
+    const PtrList<surfaceZonesInfo>& surfList,
+    const searchableSurfaces& allGeometry,
+    const labelList& surfaces
+)
+{
+    labelList unclosed(surfList.size());
+
+    label unclosedI = 0;
+    forAll(surfList, surfI)
+    {
+        if
+        (
+            surfList.set(surfI)
+         && !allGeometry[surfaces[surfI]].hasVolumeType()
+        )
+        {
+            unclosed[unclosedI++] = surfI;
+        }
+    }
+    unclosed.setSize(unclosedI);
+
+    return unclosed;
+}
+
+
 Foam::labelList Foam::surfaceZonesInfo::getAllClosedNamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList,
@@ -292,7 +315,6 @@ Foam::labelList Foam::surfaceZonesInfo::getAllClosedNamedSurfaces
 }
 
 
-// Get indices of named surfaces with a
 Foam::labelList Foam::surfaceZonesInfo::getInsidePointNamedSurfaces
 (
     const PtrList<surfaceZonesInfo>& surfList

src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H

@@ -199,6 +199,14 @@ public:
                 const labelList& surfaces
             );
 
+            //- Get indices of surfaces with a cellZone that are unclosed
+            static labelList getUnclosedNamedSurfaces
+            (
+                const PtrList<surfaceZonesInfo>& surfList,
+                const searchableSurfaces& allGeometry,
+                const labelList& surfaces
+            );
+
             //- Get indices of surfaces with a cellZone that are closed.
             static labelList getAllClosedNamedSurfaces
             (