ENH: extrudeMesh,snappyHexMesh: extrude across non-manifold boundaries

ReleaseNotes-dev

@@ -104,6 +104,11 @@
       taking film into account
     + Parallel aware
 *** *New* ptscotch decomposition method
+*** *Updated* scotch decomposition method to run in parallel by doing
+    decomposition on the master. Unfortunately scotch and ptscotch cannot
+    be linked in to the same executable.
+*** *Updated* simple decomposition method to run in parallel by doing
+    decomposition on the master.
 *** *Updated* decomposePar maps polyPatches instead of recreating them so
     polyPatches holding data can map the data.
 *** *Updated* particle tracking algorithm
@@ -183,6 +188,9 @@
       e.g. pitzDailyDirectMapped tutorial.
     + =setSet=: allows time range (e.g. 0:100) in combination with -batch argument
     to execute the commands for multiple times.
+    + =extrudeMesh=:
+            - option to add extrusion to existing mesh.
+            - works in parallel
 * Post-processing
   + =paraFoam=, =foamToVTK=: full support for polyhedral cell type in recent
      Paraview versions.

applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C

@@ -43,6 +43,7 @@ Description
 #include "addPatchCellLayer.H"
 #include "fvMesh.H"
 #include "MeshedSurfaces.H"
+#include "globalIndex.H"
 
 #include "extrudedMesh.H"
 #include "extrudeModel.H"
@@ -261,6 +262,12 @@ int main(int argc, char *argv[])
         sourceCasePath.expand();
         fileName sourceRootDir = sourceCasePath.path();
         fileName sourceCaseDir = sourceCasePath.name();
+        if (Pstream::parRun())
+        {
+            sourceCaseDir =
+                sourceCaseDir
+               /"processor" + Foam::name(Pstream::myProcNo());
+        }
         wordList sourcePatches;
         dict.lookup("sourcePatches") >> sourcePatches;
 
@@ -279,54 +286,173 @@ int main(int argc, char *argv[])
             sourceRootDir,
             sourceCaseDir
         );
+
         #include "createMesh.H"
 
         const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
 
-        // Topo change container. Either copy an existing mesh or start
-        // with empty storage (number of patches only needed for checking)
-        autoPtr<polyTopoChange> meshMod
-        (
-            (
-                mode == MESH
-              ? new polyTopoChange(mesh)
-              : new polyTopoChange(patches.size())
-            )
-        );
-
         // Extrusion engine. Either adding to existing mesh or
         // creating separate mesh.
         addPatchCellLayer layerExtrude(mesh, (mode == MESH));
 
+        const labelList meshFaces(patchFaces(patches, sourcePatches));
         indirectPrimitivePatch extrudePatch
         (
             IndirectList<face>
             (
                 mesh.faces(),
-                patchFaces(patches, sourcePatches)
+                meshFaces
             ),
             mesh.points()
         );
 
+        // Determine extrudePatch normal
+        pointField extrudePatchPointNormals
+        (
+            PatchTools::pointNormals    //calcNormals
+            (
+                mesh,
+                extrudePatch,
+                meshFaces
+            )
+        );
+
+
+        // Precalculate mesh edges for pp.edges.
+        const labelList meshEdges
+        (
+            extrudePatch.meshEdges
+            (
+                mesh.edges(),
+                mesh.pointEdges()
+            )
+        );
+
+        // Global face indices engine
+        const globalIndex globalFaces(mesh.nFaces());
+
+        // Determine extrudePatch.edgeFaces in global numbering (so across
+        // coupled patches)
+        labelListList edgeGlobalFaces
+        (
+            addPatchCellLayer::globalEdgeFaces
+            (
+                mesh,
+                globalFaces,
+                extrudePatch
+            )
+        );
+
+
+        // Determine what patches boundary edges need to get extruded into.
+        // This might actually cause edge-connected processors to become
+        // face-connected so might need to introduce new processor boundaries.
+        // Calculates:
+        //  - per pp.edge the patch to extrude into
+        //  - any additional processor boundaries (patchToNbrProc = map from
+        //    new patchID to neighbour processor)
+        //  - number of new patches (nPatches)
+
+        labelList sidePatchID;
+        label nPatches;
+        Map<label> nbrProcToPatch;
+        Map<label> patchToNbrProc;
+        addPatchCellLayer::calcSidePatch
+        (
+            mesh,
+            globalFaces,
+            edgeGlobalFaces,
+            extrudePatch,
+
+            sidePatchID,
+            nPatches,
+            nbrProcToPatch,
+            patchToNbrProc
+        );
+
+
+        // Add any patches.
+
+        label nAdded = nPatches - mesh.boundaryMesh().size();
+        reduce(nAdded, sumOp<label>());
+
+        Info<< "Adding overall " << nAdded << " processor patches." << endl;
+
+        if (nAdded > 0)
+        {
+            DynamicList<polyPatch*> newPatches(nPatches);
+            forAll(mesh.boundaryMesh(), patchI)
+            {
+                newPatches.append
+                (
+                    mesh.boundaryMesh()[patchI].clone
+                    (
+                        mesh.boundaryMesh()
+                    ).ptr()
+                );
+            }
+            for
+            (
+                label patchI = mesh.boundaryMesh().size();
+                patchI < nPatches;
+                patchI++
+            )
+            {
+                label nbrProcI = patchToNbrProc[patchI];
+                word name =
+                        "procBoundary"
+                      + Foam::name(Pstream::myProcNo())
+                      + "to"
+                      + Foam::name(nbrProcI);
+
+                Pout<< "Adding patch " << patchI
+                    << " name:" << name
+                    << " between " << Pstream::myProcNo()
+                    << " and " << nbrProcI
+                    << endl;
+
+
+                newPatches.append
+                (
+                    new processorPolyPatch
+                    (
+                        name,
+                        0,                  // size
+                        mesh.nFaces(),      // start
+                        patchI,             // index
+                        mesh.boundaryMesh(),// polyBoundaryMesh
+                        Pstream::myProcNo(),// myProcNo
+                        nbrProcI            // neighbProcNo
+                    )
+                );
+            }
+
+            // Add patches. Do no parallel updates.
+            mesh.removeFvBoundary();
+            mesh.addFvPatches(newPatches, true);
+        }
+
+
+
         // Only used for addPatchCellLayer into new mesh
-        labelList exposedPatchIDs;
+        labelList exposedPatchID;
         if (mode == PATCH)
         {
             dict.lookup("exposedPatchName") >> backPatchName;
-            exposedPatchIDs.setSize
+            exposedPatchID.setSize
             (
                 extrudePatch.size(),
                 findPatchID(patches, backPatchName)
             );
         }
 
-
+        // Determine points and extrusion
         pointField layer0Points(extrudePatch.nPoints());
         pointField displacement(extrudePatch.nPoints());
         forAll(displacement, pointI)
         {
-            const vector& patchNormal = extrudePatch.pointNormals()[pointI];
+            const vector& patchNormal = extrudePatchPointNormals[pointI];
 
             // layer0 point
             layer0Points[pointI] = model()
@@ -373,15 +499,31 @@ int main(int argc, char *argv[])
         // Layers per point
         labelList nPointLayers(extrudePatch.nPoints(), model().nLayers());
         // Displacement for first layer
-        vectorField firstLayerDisp(displacement*model().sumThickness(1));
+        vectorField firstLayerDisp = displacement*model().sumThickness(1);
+
         // Expansion ratio not used.
         scalarField ratio(extrudePatch.nPoints(), 1.0);
 
+        // Topo change container. Either copy an existing mesh or start
+        // with empty storage (number of patches only needed for checking)
+        autoPtr<polyTopoChange> meshMod
+        (
+            (
+                mode == MESH
+              ? new polyTopoChange(mesh)
+              : new polyTopoChange(patches.size())
+            )
+        );
+
         layerExtrude.setRefinement
         (
+            globalFaces,
+            edgeGlobalFaces,
+
             ratio,              // expansion ratio
             extrudePatch,       // patch faces to extrude
-            exposedPatchIDs,    // if new mesh: patches for exposed faces
+            sidePatchID,        // if boundary edge: patch to use
+            exposedPatchID,     // if new mesh: patches for exposed faces
             nFaceLayers,
             nPointLayers,
             firstLayerDisp,
@@ -401,7 +543,7 @@ int main(int argc, char *argv[])
                     model()
                     (
                         extrudePatch.localPoints()[pointI],
-                        extrudePatch.pointNormals()[pointI],
+                        extrudePatchPointNormals[pointI],
                         pPointI+1       // layer
                     )
                 );

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -43,75 +43,6 @@ defineTypeNameAndDebug(Foam::addPatchCellLayer, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Calculate global faces per pp edge.
-Foam::labelListList Foam::addPatchCellLayer::calcGlobalEdgeFaces
-(
-    const polyMesh& mesh,
-    const globalIndex& globalFaces,
-    const indirectPrimitivePatch& pp,
-    const labelList& meshEdges
-)
-{
-    //// Determine coupled edges just so we don't have to have storage
-    //// for all non-coupled edges.
-    //
-    //PackedBoolList isCoupledEdge(mesh.nEdges());
-    //
-    //const polyBoundaryMesh& patches = mesh.boundaryMesh();
-    //
-    //forAll(patches, patchI)
-    //{
-    //    const polyPatch& pp = patches[patchI];
-    //
-    //    if (pp.coupled())
-    //    {
-    //        const labelList& meshEdges = pp.meshEdges();
-    //
-    //        forAll(meshEdges, i)
-    //        {
-    //            isCoupledEdge.set(meshEdges[i], 1);
-    //        }
-    //    }
-    //}
-
-    // From mesh edge to global face labels. Sized only for pp edges.
-    labelListList globalEdgeFaces(mesh.nEdges());
-
-    const labelListList& edgeFaces = pp.edgeFaces();
-
-    forAll(edgeFaces, edgeI)
-    {
-        label meshEdgeI = meshEdges[edgeI];
-
-        //if (isCoupledEdge.get(meshEdgeI) == 1)
-        {
-            const labelList& eFaces = edgeFaces[edgeI];
-
-            // Store face and processor as unique tag.
-            labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
-            globalEFaces.setSize(eFaces.size());
-            forAll(eFaces, i)
-            {
-                globalEFaces[i] =
-                    globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
-            }
-        }
-    }
-
-    // Synchronise across coupled edges.
-    syncTools::syncEdgeList
-    (
-        mesh,
-        globalEdgeFaces,
-        uniqueEqOp(),
-        labelList()     // null value
-    );
-
-    // Extract pp part
-    return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges));
-}
-
-
 Foam::label Foam::addPatchCellLayer::nbrFace
 (
     const labelListList& edgeFaces,
@@ -316,12 +247,12 @@ Foam::labelPair Foam::addPatchCellLayer::getEdgeString
 Foam::label Foam::addPatchCellLayer::addSideFace
 (
     const indirectPrimitivePatch& pp,
-    const labelList& patchID,           // prestored patch per pp face
     const labelListList& addedCells,    // per pp face the new extruded cell
     const face& newFace,
+    const label newPatchID,
+
     const label ownFaceI,               // pp face that provides owner
     const label nbrFaceI,
-    const label patchEdgeI,             // edge to add to
     const label meshEdgeI,              // corresponding mesh edge
     const label layerI,                 // layer
     const label numEdgeFaces,           // number of layers for edge
@@ -329,8 +260,9 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     polyTopoChange& meshMod
 ) const
 {
-    // Edge to 'inflate' from
+    // Face or edge to 'inflate' from
     label inflateEdgeI = -1;
+    label inflateFaceI = -1;
 
     // Check mesh faces using edge
     if (addToMesh_)
@@ -346,8 +278,6 @@ Foam::label Foam::addPatchCellLayer::addSideFace
         }
     }
 
-    // Get my mesh face and its zone.
-    label meshFaceI = pp.addressing()[ownFaceI];
     // Zone info comes from any side patch face. Otherwise -1 since we
     // don't know what to put it in - inherit from the extruded faces?
     label zoneI = -1;   //mesh_.faceZones().whichZone(meshFaceI);
@@ -358,14 +288,15 @@ Foam::label Foam::addPatchCellLayer::addSideFace
     // Is patch edge external edge of indirectPrimitivePatch?
     if (nbrFaceI == -1)
     {
-        // External edge so external face. Patch id is obtained from
+        // External edge so external face.
-        // any other patch connected to edge.
 
         const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
         // Loop over all faces connected to edge to inflate and
-        // see if any boundary face (but not meshFaceI)
+        // see if we can find a face that is otherPatchID
-        label otherPatchID = patchID[ownFaceI];
+
+        // Get my mesh face and its zone.
+        label meshFaceI = pp.addressing()[ownFaceI];
 
         forAll(meshFaces, k)
         {
@@ -373,11 +304,14 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 
             if
             (
-                faceI != meshFaceI
+                (faceI != meshFaceI)
-             && !mesh_.isInternalFace(faceI)
+             && (patches.whichPatch(faceI) == newPatchID)
             )
             {
-                otherPatchID = patches.whichPatch(faceI);
+                // Found the patch face. Use it to inflate from
+                inflateEdgeI = -1;
+                inflateFaceI = faceI;
+
                 zoneI = mesh_.faceZones().whichZone(faceI);
                 if (zoneI != -1)
                 {
@@ -414,7 +348,7 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 
         //Pout<< "Added boundary face:" << newFace
         //    << " own:" << addedCells[ownFaceI][layerOwn]
-        //    << " patch:" << otherPatchID
+        //    << " patch:" << newPatchID
         //    << endl;
 
         addedFaceI = meshMod.setAction
@@ -426,9 +360,9 @@ Foam::label Foam::addPatchCellLayer::addSideFace
                 -1,                         // neighbour
                 -1,                         // master point
                 inflateEdgeI,               // master edge
-                -1,                         // master face
+                inflateFaceI,               // master face
                 false,                      // flux flip
-                otherPatchID,               // patch for face
+                newPatchID,                 // patch for face
                 zoneI,                      // zone for face
                 flip                        // face zone flip
             )
@@ -510,6 +444,51 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 }
 
 
+Foam::label Foam::addPatchCellLayer::findProcPatch
+(
+    const polyMesh& mesh,
+    const label nbrProcID
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(mesh.globalData().processorPatches(), i)
+    {
+        label patchI = mesh.globalData().processorPatches()[i];
+
+        if
+        (
+            refCast<const processorPolyPatch>(patches[patchI]).neighbProcNo()
+         == nbrProcID
+        )
+        {
+            return patchI;
+        }
+    }
+    return -1;
+}
+
+
+void Foam::addPatchCellLayer::setFaceProps
+(
+    const polyMesh& mesh,
+    const label faceI,
+
+    label& patchI,
+    label& zoneI,
+    bool& zoneFlip
+)
+{
+    patchI = mesh.boundaryMesh().whichPatch(faceI);
+    zoneI = mesh.faceZones().whichZone(faceI);
+    if (zoneI != -1)
+    {
+        label index = mesh.faceZones()[zoneI].whichFace(faceI);
+        zoneFlip = mesh.faceZones()[zoneI].flipMap()[index];
+    }
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from mesh
@@ -561,10 +540,251 @@ Foam::labelListList Foam::addPatchCellLayer::addedCells() const
 }
 
 
+// Calculate global faces per pp edge.
+Foam::labelListList Foam::addPatchCellLayer::globalEdgeFaces
+(
+    const polyMesh& mesh,
+    const globalIndex& globalFaces,
+    const indirectPrimitivePatch& pp
+)
+{
+    // Precalculate mesh edges for pp.edges.
+    const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
+
+    // From mesh edge to global face labels. Non-empty sublists only for
+    // pp edges.
+    labelListList globalEdgeFaces(mesh.nEdges());
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+
+    forAll(edgeFaces, edgeI)
+    {
+        label meshEdgeI = meshEdges[edgeI];
+
+        const labelList& eFaces = edgeFaces[edgeI];
+
+        // Store face and processor as unique tag.
+        labelList& globalEFaces = globalEdgeFaces[meshEdgeI];
+        globalEFaces.setSize(eFaces.size());
+        forAll(eFaces, i)
+        {
+            globalEFaces[i] = globalFaces.toGlobal(pp.addressing()[eFaces[i]]);
+        }
+    }
+
+    // Synchronise across coupled edges.
+    syncTools::syncEdgeList
+    (
+        mesh,
+        globalEdgeFaces,
+        uniqueEqOp(),
+        labelList()             // null value
+    );
+
+    // Extract pp part
+    return labelListList(UIndirectList<labelList>(globalEdgeFaces, meshEdges));
+}
+
+
+void Foam::addPatchCellLayer::calcSidePatch
+(
+    const polyMesh& mesh,
+    const globalIndex& globalFaces,
+    const labelListList& globalEdgeFaces,
+    const indirectPrimitivePatch& pp,
+
+    labelList& sidePatchID,
+    label& nPatches,
+    Map<label>& nbrProcToPatch,
+    Map<label>& patchToNbrProc
+)
+{
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    // Precalculate mesh edges for pp.edges.
+    const labelList meshEdges(pp.meshEdges(mesh.edges(), mesh.pointEdges()));
+
+    sidePatchID.setSize(pp.nEdges());
+    sidePatchID = -1;
+
+    // These also get determined but not (yet) exported:
+    // - whether face is created from other face or edge
+    // - what zone&orientation face should have
+
+    labelList inflateEdgeI(pp.nEdges(), -1);
+    labelList inflateFaceI(pp.nEdges(), -1);
+    labelList sideZoneID(pp.nEdges(), -1);
+    boolList sideFlip(pp.nEdges(), false);
+
+    nPatches = patches.size();
+
+    forAll(globalEdgeFaces, edgeI)
+    {
+        const labelList& eGlobalFaces = globalEdgeFaces[edgeI];
+        if
+        (
+            eGlobalFaces.size() == 2
+         && pp.edgeFaces()[edgeI].size() == 1
+        )
+        {
+            // Locally but not globally a boundary edge. Hence a coupled
+            // edge. Find the patch to use if on different
+            // processors.
+
+            label f0 = eGlobalFaces[0];
+            label f1 = eGlobalFaces[1];
+
+            label otherProcI = -1;
+            if (globalFaces.isLocal(f0) && !globalFaces.isLocal(f1))
+            {
+                otherProcI = globalFaces.whichProcID(f1);
+            }
+            else if (!globalFaces.isLocal(f0) && globalFaces.isLocal(f1))
+            {
+                otherProcI = globalFaces.whichProcID(f0);
+            }
+
+
+            if (otherProcI != -1)
+            {
+                sidePatchID[edgeI] = findProcPatch(mesh, otherProcI);
+                if (sidePatchID[edgeI] == -1)
+                {
+                    // Cannot find a patch to processor. See if already
+                    // marked for addition
+                    if (nbrProcToPatch.found(otherProcI))
+                    {
+                        sidePatchID[edgeI] = nbrProcToPatch[otherProcI];
+                    }
+                    else
+                    {
+                        sidePatchID[edgeI] = nPatches;
+                        nbrProcToPatch.insert(otherProcI, nPatches);
+                        patchToNbrProc.insert(nPatches, otherProcI);
+                        nPatches++;
+                    }
+                }
+            }
+        }
+    }
+
+
+
+    // Determine face properties for all other boundary edges
+    // ------------------------------------------------------
+
+    const labelListList& edgeFaces = pp.edgeFaces();
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
+        {
+            // Proper, uncoupled patch edge.
+
+            label myFaceI = pp.addressing()[edgeFaces[edgeI][0]];
+
+            // Pick up any boundary face on this edge and use its properties
+            label meshEdgeI = meshEdges[edgeI];
+            const labelList& meshFaces = mesh.edgeFaces()[meshEdgeI];
+
+            forAll(meshFaces, k)
+            {
+                label faceI = meshFaces[k];
+
+                if (faceI != myFaceI && !mesh.isInternalFace(faceI))
+                {
+                    setFaceProps
+                    (
+                        mesh,
+                        faceI,
+
+                        sidePatchID[edgeI],
+                        sideZoneID[edgeI],
+                        sideFlip[edgeI]
+                    );
+                    inflateFaceI[edgeI] = faceI;
+                    inflateEdgeI[edgeI] = -1;
+
+                    break;
+                }
+            }
+        }
+    }
+
+
+
+    // Now hopefully every boundary edge has a side patch. Check
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && sidePatchID[edgeI] == -1)
+        {
+            const edge& e = pp.edges()[edgeI];
+            FatalErrorIn("addPatchCellLayer::calcSidePatch(..)")
+                << "Have no sidePatchID for edge " << edgeI << " points "
+                << pp.points()[pp.meshPoints()[e[0]]]
+                << pp.points()[pp.meshPoints()[e[1]]]
+                << abort(FatalError);
+        }
+    }
+
+
+
+    // Now we have sidepatch see if we have patchface or edge to inflate
+    // from.
+    forAll(edgeFaces, edgeI)
+    {
+        if (edgeFaces[edgeI].size() == 1 && inflateFaceI[edgeI] == -1)
+        {
+            // 1. Do we have a boundary face to inflate from
+
+            label myFaceI = pp.addressing()[edgeFaces[edgeI][0]];
+
+            // Pick up any boundary face on this edge and use its properties
+            label meshEdgeI = meshEdges[edgeI];
+            const labelList& meshFaces = mesh.edgeFaces()[meshEdgeI];
+
+            forAll(meshFaces, k)
+            {
+                label faceI = meshFaces[k];
+
+                if (faceI != myFaceI)
+                {
+                    if (mesh.isInternalFace(faceI))
+                    {
+                        inflateEdgeI[edgeI] = meshEdgeI;
+                    }
+                    else
+                    {
+                        if (patches.whichPatch(faceI) == sidePatchID[edgeI])
+                        {
+                            setFaceProps
+                            (
+                                mesh,
+                                faceI,
+
+                                sidePatchID[edgeI],
+                                sideZoneID[edgeI],
+                                sideFlip[edgeI]
+                            );
+                            inflateFaceI[edgeI] = faceI;
+                            inflateEdgeI[edgeI] = -1;
+
+                            break;
+                        }
+                    }
+                }
+            }
+        }
+    }
+}
+
+
 void Foam::addPatchCellLayer::setRefinement
 (
+    const globalIndex& globalFaces,
+    const labelListList& globalEdgeFaces,
     const scalarField& expansionRatio,
     const indirectPrimitivePatch& pp,
+    const labelList& sidePatchID,
     const labelList& exposedPatchID,
     const labelList& nFaceLayers,
     const labelList& nPointLayers,
@@ -575,7 +795,7 @@ void Foam::addPatchCellLayer::setRefinement
     if (debug)
     {
         Pout<< "addPatchCellLayer::setRefinement : Adding up to "
-            << max(nPointLayers)
+            << gMax(nPointLayers)
             << " layers of cells to indirectPrimitivePatch with "
             << pp.nPoints() << " points" << endl;
     }
@@ -788,8 +1008,6 @@ void Foam::addPatchCellLayer::setRefinement
 
                 label meshEdgeI = meshEdges[edgeI];
 
-                // Mesh faces using edge
-
                 // Mesh faces using edge
                 const labelList& meshFaces = mesh_.edgeFaces(meshEdgeI, ef);
 
@@ -1213,22 +1431,6 @@ void Foam::addPatchCellLayer::setRefinement
     }
 
 
-    // Global indices engine
-    const globalIndex globalFaces(mesh_.nFaces());
-
-    // Get for all pp edgeFaces a unique faceID
-    labelListList globalEdgeFaces
-    (
-         calcGlobalEdgeFaces
-         (
-            mesh_,
-            globalFaces,
-            pp,
-            meshEdges
-        )
-    );
-
-
     // Mark off which edges have been extruded
     boolList doneEdge(pp.nEdges(), false);
 
@@ -1474,16 +1676,17 @@ void Foam::addPatchCellLayer::setRefinement
                         addSideFace
                         (
                             pp,
-                            patchID,
                             addedCells,
-                            newFace,
+
+                            newFace,                // vertices of new face
+                            sidePatchID[startEdgeI],// -1 or patch for face
+
                             patchFaceI,
                             nbrFaceI,
-                            startEdgeI,     // edge to inflate from
+                            meshEdgeI,          // (mesh) edge to inflate
-                            meshEdgeI,      // corresponding mesh edge
+                            i,                  // layer
-                            i,
+                            numEdgeSideFaces,   // num layers
-                            numEdgeSideFaces,
+                            meshFaces,          // edgeFaces
-                            meshFaces,
                             meshMod
                         );
                     }

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.H

@@ -181,16 +181,6 @@ class addPatchCellLayer
 
     // Private Member Functions
 
-        //- Per patch edge the pp faces (in global indices) using it. Uses
-        //  uniqueEqOp() to remove duplicates.
-        labelListList calcGlobalEdgeFaces
-        (
-            const polyMesh& mesh,
-            const globalIndex& globalFaces,
-            const indirectPrimitivePatch& pp,
-            const labelList& meshEdges
-        );
-
         //- Get the face on the other side of the edge.
         static label nbrFace
         (
@@ -226,12 +216,13 @@ class addPatchCellLayer
         label addSideFace
         (
             const indirectPrimitivePatch&,
-            const labelList& patchID,
             const labelListList& addedCells,
+
             const face& newFace,
+            const label newPatchID,
+
             const label ownFaceI,
             const label nbrFaceI,
-            const label patchEdgeI,
             const label meshEdgeI,
             const label layerI,
             const label numEdgeFaces,
@@ -239,6 +230,18 @@ class addPatchCellLayer
             polyTopoChange&
         ) const;
 
+        //- Find patch to neighbouring processor
+        static label findProcPatch(const polyMesh&, const label nbrProcID);
+
+        //- Extract properties from mesh face
+        static void setFaceProps
+        (
+            const polyMesh&,
+            const label,
+            label&,
+            label&,
+            bool&
+        );
 
         //- Disallow default bitwise copy construct
         addPatchCellLayer(const addPatchCellLayer&);
@@ -256,7 +259,7 @@ public:
     // Constructors
 
         //- Construct from mesh.
-        addPatchCellLayer(const polyMesh& mesh, const bool addToMesh = true);
+        addPatchCellLayer(const polyMesh&, const bool addToMesh = true);
 
 
     // Member Functions
@@ -291,6 +294,33 @@ public:
 
         // Edit
 
+            //- Per patch edge the pp faces (in global indices) using it. Uses
+            //  uniqueEqOp() to remove duplicates.
+            static labelListList globalEdgeFaces
+            (
+                const polyMesh&,
+                const globalIndex& globalFaces,
+                const indirectPrimitivePatch& pp
+            );
+
+            //- Boundary edges get extruded into boundary faces. Determine patch
+            //  for these faces. This might be a to-be-created processor patch
+            //  (patchI >= mesh.boundaryMesh().size()) in which case the
+            //  nbrProcToPatch, patchToNbrProc give the correspondence. nPatches
+            //  is the new number of patches.
+            static void calcSidePatch
+            (
+                const polyMesh&,
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
+                const indirectPrimitivePatch& pp,
+
+                labelList& sidePatchID,
+                label& nPatches,
+                Map<label>& nbrProcToPatch,
+                Map<label>& patchToNbrProc
+            );
+
             //- Play commands into polyTopoChange to create layers on top
             //  of indirectPrimitivePatch (have to be outside faces).
             //  Gets displacement per patch point.
@@ -313,8 +343,11 @@ public:
             //        (instead of e.g. from patch faces)
             void setRefinement
             (
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
                 const scalarField& expansionRatio,
                 const indirectPrimitivePatch& pp,
+                const labelList& sidePatchID,
                 const labelList& exposedPatchID,
                 const labelList& nFaceLayers,
                 const labelList& nPointLayers,
@@ -326,20 +359,26 @@ public:
             //- Add with constant expansion ratio and same nLayers everywhere
             void setRefinement
             (
+                const globalIndex& globalFaces,
+                const labelListList& globalEdgeFaces,
                 const label nLayers,
                 const indirectPrimitivePatch& pp,
+                const labelList& sidePatchID,
                 const vectorField& overallDisplacement,
                 polyTopoChange& meshMod
             )
             {
                 setRefinement
                 (
-                    scalarField(pp.nPoints(), 1.0), // expansion ration
+                    globalFaces,
+                    globalEdgeFaces,
+                    scalarField(pp.nPoints(), 1.0),     // expansion ration
                     pp,
+                    sidePatchID,
                     labelList(0),
-                    labelList(pp.size(), nLayers),
+                    labelList(pp.size(), nLayers),      // nFaceLayers
-                    labelList(pp.nPoints(), nLayers),
+                    labelList(pp.nPoints(), nLayers),   // nPointLayers
-                    overallDisplacement / nLayers,
+                    overallDisplacement / nLayers,      // firstLayerDisp
                     meshMod
                 );
             }

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ Description
 #include "layerParameters.H"
 #include "combineFaces.H"
 #include "IOmanip.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -888,68 +889,6 @@ void Foam::autoLayerDriver::handleNonManifolds
     Info<< "Outside of local patch is multiply connected across edges or"
         << " points at " << nNonManif << " points." << endl;
 
-
-    const labelList& meshPoints = pp.meshPoints();
-
-
-    // 2. Parallel check
-    // For now disable extrusion at any shared edges.
-    const labelHashSet sharedEdgeSet(mesh.globalData().sharedEdgeLabels());
-
-    forAll(pp.edges(), edgeI)
-    {
-        if (sharedEdgeSet.found(meshEdges[edgeI]))
-        {
-            const edge& e = mesh.edges()[meshEdges[edgeI]];
-
-            Pout<< "Disabling extrusion at edge "
-                << mesh.points()[e[0]] << mesh.points()[e[1]]
-                << " since it is non-manifold across coupled patches."
-                << endl;
-
-            nonManifoldPoints.insert(e[0]);
-            nonManifoldPoints.insert(e[1]);
-        }
-    }
-
-    // 3b. extrusion can produce multiple faces between 2 cells
-    // across processor boundary
-    // This occurs when a coupled face shares more than 1 edge with a
-    // non-coupled boundary face.
-    // This is now correctly handled by addPatchCellLayer in that it
-    // extrudes a single face from the stringed up edges.
-
-
-    nNonManif = returnReduce(nonManifoldPoints.size(), sumOp<label>());
-
-    if (nNonManif > 0)
-    {
-        Info<< "Outside of patches is multiply connected across edges or"
-            << " points at " << nNonManif << " points." << nl
-            << "Writing " << nNonManif
-            << " points where this happens to pointSet "
-            << nonManifoldPoints.name() << nl
-            << "and setting layer thickness to zero on these points."
-            << endl;
-
-        nonManifoldPoints.instance() = mesh.time().timeName();
-        nonManifoldPoints.write();
-
-        forAll(meshPoints, patchPointI)
-        {
-            if (nonManifoldPoints.found(meshPoints[patchPointI]))
-            {
-                unmarkExtrusion
-                (
-                    patchPointI,
-                    patchDisp,
-                    patchNLayers,
-                    extrudeStatus
-                );
-            }
-        }
-    }
-
     Info<< "Set displacement to zero for all " << nNonManif
         << " non-manifold points" << endl;
 }
@@ -1348,7 +1287,6 @@ void Foam::autoLayerDriver::setNumLayers
 }
 
 
-// Grow no-extrusion layer
 void Foam::autoLayerDriver::growNoExtrusion
 (
     const indirectPrimitivePatch& pp,
@@ -1439,6 +1377,103 @@ void Foam::autoLayerDriver::growNoExtrusion
 }
 
 
+void Foam::autoLayerDriver::determineSidePatches
+(
+    const globalIndex& globalFaces,
+    const labelListList& edgeGlobalFaces,
+    const indirectPrimitivePatch& pp,
+
+    labelList& sidePatchID
+)
+{
+    // Sometimes edges-to-be-extruded are on more than 2 processors.
+    // Work out which 2 hold the faces to be extruded and thus which procpatch
+    // the side-face should be in. As an additional complication this might
+    // mean that 2 procesors that were only edge-connected now suddenly need
+    // to become face-connected i.e. have a processor patch between them.
+
+    fvMesh& mesh = meshRefiner_.mesh();
+
+    // Determine sidePatchID. Any additional processor boundary gets added to
+    // patchToNbrProc,nbrProcToPatch and nPatches gets set to the new number
+    // of patches.
+    label nPatches;
+    Map<label> nbrProcToPatch;
+    Map<label> patchToNbrProc;
+    addPatchCellLayer::calcSidePatch
+    (
+        mesh,
+        globalFaces,
+        edgeGlobalFaces,
+        pp,
+
+        sidePatchID,
+        nPatches,
+        nbrProcToPatch,
+        patchToNbrProc
+    );
+
+    label nOldPatches = mesh.boundaryMesh().size();
+    label nAdded = returnReduce(nPatches-nOldPatches, sumOp<label>());
+    Info<< nl << "Adding in total " << nAdded/2 << " inter-processor patches to"
+        << " handle extrusion of non-manifold processor boundaries."
+        << endl;
+
+    if (nAdded > 0)
+    {
+        // We might not add patches in same order as in patchToNbrProc
+        // so prepare to renumber sidePatchID
+        Map<label> wantedToAddedPatch;
+
+        for (label patchI = nOldPatches; patchI < nPatches; patchI++)
+        {
+            label nbrProcI = patchToNbrProc[patchI];
+            word name =
+                    "procBoundary"
+                  + Foam::name(Pstream::myProcNo())
+                  + "to"
+                  + Foam::name(nbrProcI);
+
+            dictionary patchDict;
+            patchDict.add("type", processorPolyPatch::typeName);
+            patchDict.add("myProcNo", Pstream::myProcNo());
+            patchDict.add("neighbProcNo", nbrProcI);
+            patchDict.add("nFaces", 0);
+            patchDict.add("startFace", mesh.nFaces());
+
+            Pout<< "Adding patch " << patchI
+                << " name:" << name
+                << " between " << Pstream::myProcNo()
+                << " and " << nbrProcI
+                << endl;
+
+            label procPatchI = meshRefiner_.appendPatch
+            (
+                mesh,
+                mesh.boundaryMesh().size(), // new patch index
+                name,
+                patchDict
+            );
+            wantedToAddedPatch.insert(patchI, procPatchI);
+        }
+
+        // Renumber sidePatchID
+        forAll(sidePatchID, i)
+        {
+            label patchI = sidePatchID[i];
+            Map<label>::const_iterator fnd = wantedToAddedPatch.find(patchI);
+            if (fnd != wantedToAddedPatch.end())
+            {
+                sidePatchID[i] = fnd();
+            }
+        }
+
+        mesh.clearOut();
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh()).updateMesh();
+    }
+}
+
+
 void Foam::autoLayerDriver::calculateLayerThickness
 (
     const indirectPrimitivePatch& pp,
@@ -2553,6 +2588,37 @@ void Foam::autoLayerDriver::addLayers
         )
     );
 
+
+    // Global face indices engine
+    const globalIndex globalFaces(mesh.nFaces());
+
+    // Determine extrudePatch.edgeFaces in global numbering (so across
+    // coupled patches). This is used only to string up edges between coupled
+    // faces (all edges between same (global)face indices get extruded).
+    labelListList edgeGlobalFaces
+    (
+        addPatchCellLayer::globalEdgeFaces
+        (
+            mesh,
+            globalFaces,
+            pp
+        )
+    );
+
+    // Determine patches for extruded boundary edges. Calculates if any
+    // additional processor patches need to be constructed.
+
+    labelList sidePatchID;
+    determineSidePatches
+    (
+        globalFaces,
+        edgeGlobalFaces,
+        pp,
+
+        sidePatchID
+    );
+
+
     // Construct iterative mesh mover.
     Info<< "Constructing mesh displacer ..." << endl;
 
@@ -3038,8 +3104,12 @@ void Foam::autoLayerDriver::addLayers
         // Not add layer if patchDisp is zero.
         addLayer.setRefinement
         (
+            globalFaces,
+            edgeGlobalFaces,
+
             invExpansionRatio,
             pp(),
+            sidePatchID,        // boundary patch for extruded boundary edges
             labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -243,6 +243,16 @@ class autoLayerDriver
                     List<extrudeMode>& extrudeStatus
                 ) const;
 
+                //- See what patches boundaryedges should be extruded into
+                void determineSidePatches
+                (
+                    const globalIndex& globalFaces,
+                    const labelListList& edgeGlobalFaces,
+                    const indirectPrimitivePatch& pp,
+
+                    labelList& sidePatchID
+                );
+
                 //- Calculate pointwise wanted and minimum thickness.
                 //  thickness: wanted thickness
                 //  minthickness: when to give up and not extrude

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -1549,60 +1549,28 @@ void Foam::meshRefinement::checkCoupledFaceZones(const polyMesh& mesh)
 }
 
 
-// Adds patch if not yet there. Returns patchID.
+Foam::label Foam::meshRefinement::appendPatch
-Foam::label Foam::meshRefinement::addPatch
 (
     fvMesh& mesh,
+    const label insertPatchI,
     const word& patchName,
-    const dictionary& patchInfo
+    const dictionary& patchDict
 )
 {
-    polyBoundaryMesh& polyPatches =
-        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
-
-    const label patchI = polyPatches.findPatchID(patchName);
-    if (patchI != -1)
-    {
-        // Already there
-        return patchI;
-    }
-
-
-    label insertPatchI = polyPatches.size();
-    label startFaceI = mesh.nFaces();
-
-    forAll(polyPatches, patchI)
-    {
-        const polyPatch& pp = polyPatches[patchI];
-
-        if (isA<processorPolyPatch>(pp))
-        {
-            insertPatchI = patchI;
-            startFaceI = pp.start();
-            break;
-        }
-    }
-
-
-    // Below is all quite a hack. Feel free to change once there is a better
-    // mechanism to insert and reorder patches.
-
     // Clear local fields and e.g. polyMesh parallelInfo.
     mesh.clearOut();
 
-    label sz = polyPatches.size();
+    polyBoundaryMesh& polyPatches =
-
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
     fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
 
-    dictionary patchDict(patchInfo);
+    label patchI = polyPatches.size();
-    patchDict.set("nFaces", 0);
-    patchDict.set("startFace", startFaceI);
 
     // Add polyPatch at the end
-    polyPatches.setSize(sz+1);
+    polyPatches.setSize(patchI+1);
     polyPatches.set
     (
-        sz,
+        patchI,
         polyPatch::New
         (
             patchName,
@@ -1611,13 +1579,13 @@ Foam::label Foam::meshRefinement::addPatch
             polyPatches
         )
     );
-    fvPatches.setSize(sz+1);
+    fvPatches.setSize(patchI+1);
     fvPatches.set
     (
-        sz,
+        patchI,
         fvPatch::New
         (
-            polyPatches[sz],  // point to newly added polyPatch
+            polyPatches[patchI],  // point to newly added polyPatch
             mesh.boundary()
         )
     );
@@ -1675,21 +1643,69 @@ Foam::label Foam::meshRefinement::addPatch
         mesh,
         calculatedFvPatchField<tensor>::typeName
     );
+    return patchI;
+}
+
+
+// Adds patch if not yet there. Returns patchID.
+Foam::label Foam::meshRefinement::addPatch
+(
+    fvMesh& mesh,
+    const word& patchName,
+    const dictionary& patchInfo
+)
+{
+    polyBoundaryMesh& polyPatches =
+        const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
+    fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
+
+    const label patchI = polyPatches.findPatchID(patchName);
+    if (patchI != -1)
+    {
+        // Already there
+        return patchI;
+    }
+
+
+    label insertPatchI = polyPatches.size();
+    label startFaceI = mesh.nFaces();
+
+    forAll(polyPatches, patchI)
+    {
+        const polyPatch& pp = polyPatches[patchI];
+
+        if (isA<processorPolyPatch>(pp))
+        {
+            insertPatchI = patchI;
+            startFaceI = pp.start();
+            break;
+        }
+    }
+
+    dictionary patchDict(patchInfo);
+    patchDict.set("nFaces", 0);
+    patchDict.set("startFace", startFaceI);
+
+    // Below is all quite a hack. Feel free to change once there is a better
+    // mechanism to insert and reorder patches.
+
+    label addedPatchI = appendPatch(mesh, insertPatchI, patchName, patchDict);
+
 
     // Create reordering list
     // patches before insert position stay as is
-    labelList oldToNew(sz+1);
+    labelList oldToNew(addedPatchI+1);
     for (label i = 0; i < insertPatchI; i++)
     {
         oldToNew[i] = i;
     }
     // patches after insert position move one up
-    for (label i = insertPatchI; i < sz; i++)
+    for (label i = insertPatchI; i < addedPatchI; i++)
     {
         oldToNew[i] = i+1;
     }
     // appended patch gets moved to insert position
-    oldToNew[sz] = insertPatchI;
+    oldToNew[addedPatchI] = insertPatchI;
 
     // Shuffle into place
     polyPatches.reorder(oldToNew);

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -747,8 +747,17 @@ public:
 
         // Other topo changes
 
+            //- Helper:append patch to end of mesh.
+            static label appendPatch
+            (
+                fvMesh&,
+                const label insertPatchI,
+                const word&,
+                const dictionary&
+            );
+
             //- Helper:add patch to mesh. Update all registered fields.
-            //  Use addMeshedPatch to add patches originating from surfaces.
+            //  Used by addMeshedPatch to add patches originating from surfaces.
             static label addPatch(fvMesh&, const word& name, const dictionary&);
 
             //- Add patch originating from meshing. Update meshedPatches_.

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -1007,6 +1007,82 @@ void Foam::refinementSurfaces::findNearestRegion
 }
 
 
+void Foam::refinementSurfaces::findNearestRegion
+(
+    const labelList& surfacesToTest,
+    const pointField& samples,
+    const scalarField& nearestDistSqr,
+    labelList& hitSurface,
+    List<pointIndexHit>& hitInfo,
+    labelList& hitRegion,
+    vectorField& hitNormal
+) const
+{
+    labelList geometries(UIndirectList<label>(surfaces_, surfacesToTest));
+
+    // Do the tests. Note that findNearest returns index in geometries.
+    searchableSurfacesQueries::findNearest
+    (
+        allGeometry_,
+        geometries,
+        samples,
+        nearestDistSqr,
+        hitSurface,
+        hitInfo
+    );
+
+    // Rework the hitSurface to be surface (i.e. index into surfaces_)
+    forAll(hitSurface, pointI)
+    {
+        if (hitSurface[pointI] != -1)
+        {
+            hitSurface[pointI] = surfacesToTest[hitSurface[pointI]];
+        }
+    }
+
+    // Collect the region
+    hitRegion.setSize(hitSurface.size());
+    hitRegion = -1;
+    hitNormal.setSize(hitSurface.size());
+    hitNormal = vector::zero;
+
+    forAll(surfacesToTest, i)
+    {
+        label surfI = surfacesToTest[i];
+
+        // Collect hits for surfI
+        const labelList localIndices(findIndices(hitSurface, surfI));
+
+        List<pointIndexHit> localHits
+        (
+            UIndirectList<pointIndexHit>
+            (
+                hitInfo,
+                localIndices
+            )
+        );
+
+        // Region
+        labelList localRegion;
+        allGeometry_[surfaces_[surfI]].getRegion(localHits, localRegion);
+
+        forAll(localIndices, i)
+        {
+            hitRegion[localIndices[i]] = localRegion[i];
+        }
+
+        // Normal
+        vectorField localNormal;
+        allGeometry_[surfaces_[surfI]].getNormal(localHits, localNormal);
+
+        forAll(localIndices, i)
+        {
+            hitNormal[localIndices[i]] = localNormal[i];
+        }
+    }
+}
+
+
 //// Find intersection with max of edge. Return -1 or the surface
 //// with the highest maxLevel above currentLevel
 //Foam::label Foam::refinementSurfaces::findHighestIntersection

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -327,6 +327,19 @@ public:
                 labelList& hitRegion
             ) const;
 
+            //- Find nearest point on surfaces. Return surface, region and
+            //  normal on surface (so not global surface)
+            void findNearestRegion
+            (
+                const labelList& surfacesToTest,
+                const pointField& samples,
+                const scalarField& nearestDistSqr,
+                labelList& hitSurface,
+                List<pointIndexHit>& hitInfo,
+                labelList& hitRegion,
+                vectorField& hitNormal
+            ) const;
+
             //- Detect if a point is 'inside' (closed) surfaces.
             //  Returns -1 if not, returns first surface it is.
             void findInside