BUG: Preserve face zone orientation (fvMeshSubset, removeCels, addPatchCellLayer)

src/dynamicMesh/polyTopoChange/polyTopoChange/addPatchCellLayer.C

@@ -349,7 +349,10 @@ Foam::label Foam::addPatchCellLayer::addSideFace
 
     // Get my mesh face and its zone.
     label meshFaceI = pp.addressing()[ownFaceI];
-    label zoneI = mesh_.faceZones().whichZone(meshFaceI);
+    // 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);
+    bool flip = false;
 
     label addedFaceI = -1;
 
@@ -376,6 +379,12 @@ Foam::label Foam::addPatchCellLayer::addSideFace
             )
             {
                 otherPatchID = patches.whichPatch(faceI);
+                zoneI = mesh_.faceZones().whichZone(faceI);
+                if (zoneI != -1)
+                {
+                    label index = mesh_.faceZones()[zoneI].whichFace(faceI);
+                    flip = mesh_.faceZones()[zoneI].flipMap()[index];
+                }
                 break;
             }
         }
@@ -422,7 +431,7 @@ Foam::label Foam::addPatchCellLayer::addSideFace
                 false,                      // flux flip
                 otherPatchID,               // patch for face
                 zoneI,                      // zone for face
-                false                       // face zone flip
+                flip                        // face zone flip
             )
         );
     }
@@ -488,7 +497,7 @@ Foam::label Foam::addPatchCellLayer::addSideFace
                 false,                      // flux flip
                 -1,                         // patch for face
                 zoneI,                      // zone for face
-                false                       // face zone flip
+                flip                        // face zone flip
             )
         );
 
@@ -1027,12 +1036,21 @@ void Foam::addPatchCellLayer::setRefinement
                 // Get new neighbour
                 label nei;
                 label patchI;
+                label zoneI = -1;
+                bool flip = false;
+
 
                 if (i == addedCells[patchFaceI].size()-1)
                 {
                     // Top layer so is patch face.
                     nei = -1;
                     patchI = patchID[patchFaceI];
+                    zoneI = mesh_.faceZones().whichZone(meshFaceI);
+                    if (zoneI != -1)
+                    {
+                        const faceZone& fz = mesh_.faceZones()[zoneI];
+                        flip = fz.flipMap()[fz.whichFace(meshFaceI)];
+                    }
                 }
                 else
                 {
@@ -1055,7 +1073,7 @@ void Foam::addPatchCellLayer::setRefinement
                         false,                      // flux flip
                         patchI,                     // patch for face
                         zoneI,                      // zone for face
-                        false                       // face zone flip
+                        flip                        // face zone flip
                     )
                 );
             }
@@ -1076,8 +1094,6 @@ void Foam::addPatchCellLayer::setRefinement
 
                 layerFaces_[patchFaceI][0] = meshFaceI;
 
-                label zoneI = mesh_.faceZones().whichZone(meshFaceI);
-
                 meshMod.setAction
                 (
                     polyModifyFace
@@ -1088,8 +1104,8 @@ void Foam::addPatchCellLayer::setRefinement
                         addedCells[patchFaceI][0],      // neighbour
                         false,                          // face flip
                         -1,                             // patch for face
-                        false,                          // remove from zone
+                        true, //false,                  // remove from zone
-                        zoneI,                          // zone for face
+                        -1, //zoneI,                    // zone for face
                         false                           // face flip in zone
                     )
                 );

src/dynamicMesh/polyTopoChange/polyTopoChange/removeCells.C

@@ -322,7 +322,9 @@ void Foam::removeCells::setRefinement
                 if (zoneID >= 0)
                 {
                     const faceZone& fZone = faceZones[zoneID];
-                    zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
+                    // Note: we reverse the owner/neighbour of the face
+                    // so should also select the other side of the zone
+                    zoneFlip = !fZone.flipMap()[fZone.whichFace(faceI)];
                 }
 
                 //Pout<< "Putting exposed internal face " << faceI

src/finiteVolume/fvMesh/fvMeshSubset/fvMeshSubset.C

@@ -276,28 +276,50 @@ void Foam::fvMeshSubset::subsetZones()
     {
         const faceZone& fz = faceZones[i];
 
-        // Create list of mesh faces part of the new zone
+        // Expand faceZone to full mesh
-        labelList subAddressing
+        // +1 : part of faceZone, flipped
-        (
+        // -1 :    ,,           , unflipped
-            subset
+        //  0 : not part of faceZone
-            (
+        labelList zone(baseMesh().nFaces(), 0);
-                baseMesh().nFaces(),
-                fz,
-                faceMap()
-            )
-        );
-
-        // Flipmap for all mesh faces
-        boolList fullFlipStatus(baseMesh().nFaces(), false);
         forAll(fz, j)
         {
-            fullFlipStatus[fz[j]] = fz.flipMap()[j];
+            if (fz.flipMap()[j])
-        }
-        // Extract sub part
-        boolList subFlipStatus(subAddressing.size(), false);
-        forAll(subAddressing, j)
             {
-            subFlipStatus[j] = fullFlipStatus[faceMap()[subAddressing[j]]];
+                zone[fz[j]] = 1;
+            }
+            else
+            {
+                zone[fz[j]] = -1;
+            }
+        }
+
+        // Select faces
+        label nSub = 0;
+        forAll(faceMap(), j)
+        {
+            if (zone[faceMap()[j]] != 0)
+            {
+                nSub++;
+            }
+        }
+        labelList subAddressing(nSub);
+        boolList subFlipStatus(nSub);
+        nSub = 0;
+        forAll(faceMap(), subFaceI)
+        {
+            label meshFaceI = faceMap()[subFaceI];
+            if (zone[meshFaceI] != 0)
+            {
+                subAddressing[nSub] = subFaceI;
+                label subOwner = subMesh().faceOwner()[subFaceI];
+                label baseOwner = baseMesh().faceOwner()[meshFaceI];
+                // If subowner is the same cell as the base keep the flip status
+                bool sameOwner = (cellMap()[subOwner] == baseOwner);
+                bool flip = (zone[meshFaceI] == 1);
+                subFlipStatus[nSub] = (sameOwner == flip);
+
+                nSub++;
+            }
         }
 
         fZonePtrs[i] = new faceZone