meshRefinementBaffles: Correct faceZone orientation

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1479
Patches provided by Bruno Santos based on the work of Mattijs Janssens

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

@@ -369,7 +369,7 @@ private:
                 const labelList& globalToSlavePatch
             ) const;
 
-            //- Determine patches for baffles
+            //- Determine patches for baffles on all intersected unnamed faces
             void getBafflePatches
             (
                 const labelList& globalToMasterPatch,
@@ -514,6 +514,8 @@ private:
                 labelList& cellToZone
             ) const;
 
+            //- Make namedSurfaceIndex consistent with cellToZone
+            //  - clear out any blocked faces inbetween same cell zone.
             void makeConsistentFaceIndex
             (
                 const labelList& cellToZone,
@@ -536,11 +538,11 @@ private:
 
             // Some patch utilities
 
-            //- Get all faces in namedSurfaceIndex that have no cellZone on
+            //- Get all faces in faceToZone that have no cellZone on
             //  either side.
             labelList freeStandingBaffleFaces
             (
-                const labelList& namedSurfaceIndex,
+                const labelList& faceToZone,
                 const labelList& cellToZone,
                 const labelList& neiCellZone
             ) const;

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

@@ -51,9 +51,6 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Repatches external face or creates baffle for internal face
-// with user specified patches (might be different for both sides).
-// Returns label of added face.
 Foam::label Foam::meshRefinement::createBaffle
 (
     const label faceI,
@@ -128,52 +125,11 @@ Foam::label Foam::meshRefinement::createBaffle
 }
 
 
-//// Check if we are a boundary face and normal of surface does
-//// not align with test vector. In this case there'd probably be
-//// a freestanding 'baffle' so we might as well not create it.
-//// Note that since it is not a proper baffle we cannot detect it
-//// afterwards so this code cannot be merged with the
-//// filterDuplicateFaces code.
-//bool Foam::meshRefinement::validBaffleTopology
-//(
-//    const label faceI,
-//    const vector& n1,
-//    const vector& n2,
-//    const vector& testDir
-//) const
-//{
-//
-//    label patchI = mesh_.boundaryMesh().whichPatch(faceI);
-//    if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
-//    {
-//        return true;
-//    }
-//    else if (mag(n1&n2) > cos(degToRad(30)))
-//    {
-//        // Both normals aligned. Check that test vector perpendicularish to
-//        // surface normal
-//        scalar magTestDir = mag(testDir);
-//        if (magTestDir > VSMALL)
-//        {
-//            if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45)))
-//            {
-//                //Pout<< "** disabling baffling face "
-//                //    << mesh_.faceCentres()[faceI] << endl;
-//                return false;
-//            }
-//        }
-//    }
-//    return true;
-//}
-
-
-// Determine patches for baffles on all intersected unnamed faces
 void Foam::meshRefinement::getBafflePatches
 (
     const labelList& globalToMasterPatch,
     const labelList& neiLevel,
     const pointField& neiCc,
-
     labelList& ownPatch,
     labelList& neiPatch
 ) const
@@ -386,7 +342,6 @@ Foam::Map<Foam::labelPair> Foam::meshRefinement::getZoneBafflePatches
 }
 
 
-// Create baffle for every face where ownPatch != -1
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
 (
     const labelList& ownPatch,
@@ -659,17 +614,16 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
 }
 
 
-// Extract those baffles (duplicate) faces that are on the edge of a baffle
-// region. These are candidates for merging.
-// Done by counting the number of baffles faces per mesh edge. If edge
-// has 2 boundary faces and both are baffle faces it is the edge of a baffle
-// region.
 Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
 (
     const List<labelPair>& couples,
     const scalar planarAngle
 ) const
 {
+    // Done by counting the number of baffles faces per mesh edge. If edge
+    // has 2 boundary faces and both are baffle faces it is the edge of a baffle
+    // region.
+
     // All duplicate faces on edge of the patch are to be merged.
     // So we count for all edges of duplicate faces how many duplicate
     // faces use them.
@@ -918,7 +872,6 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
 }
 
 
-// Merge baffles. Gets pairs of faces.
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
 (
     const List<labelPair>& couples
@@ -1420,10 +1373,6 @@ bool Foam::meshRefinement::calcRegionToZone
 }
 
 
-// Finds region per cell. Assumes:
-// - region containing keepPoint does not go into a cellZone
-// - all other regions can be found by crossing faces marked in
-//   namedSurfaceIndex.
 void Foam::meshRefinement::findCellZoneTopo
 (
     const point& keepPoint,
@@ -1432,6 +1381,11 @@ void Foam::meshRefinement::findCellZoneTopo
     labelList& cellToZone
 ) const
 {
+    // Assumes:
+    // - region containing keepPoint does not go into a cellZone
+    // - all other regions can be found by crossing faces marked in
+    //   namedSurfaceIndex.
+
     // Analyse regions. Reuse regionsplit
     boolList blockedFace(mesh_.nFaces());
 
@@ -1629,8 +1583,6 @@ void Foam::meshRefinement::findCellZoneTopo
 }
 
 
-// Make namedSurfaceIndex consistent with cellToZone - clear out any blocked
-// faces inbetween same cell zone.
 void Foam::meshRefinement::makeConsistentFaceIndex
 (
     const labelList& cellToZone,
@@ -1810,7 +1762,7 @@ void Foam::meshRefinement::handleSnapProblems
 
 Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
 (
-    const labelList& namedSurfaceIndex,
+    const labelList& faceToZone,
     const labelList& cellToZone,
     const labelList& neiCellZone
 ) const
@@ -1820,17 +1772,30 @@ Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
     const labelList& faceNeighbour = mesh_.faceNeighbour();
 
 
-    DynamicList<label> faceLabels(mesh_.nFaces()/20);
+    // We want to pick up the faces to orient. These faces come in
+    // two variants:
+    // - faces originating from stand-alone faceZones
+    //   (these will most likely have no cellZone on either side so
+    //    ownZone and neiZone both -1)
+    // - sticky-up faces originating from a 'bulge' in a outside of
+    //   a cellZone. These will have the same cellZone on either side.
+    //   How to orient these is not really clearly defined so do them
+    //   same as stand-alone faceZone faces for now. (Normally these will
+    //   already have been removed by the 'allowFreeStandingZoneFaces=false'
+    //   default setting)
+
+    // Note that argument neiCellZone will have -1 on uncoupled boundaries.
+
+    DynamicList<label> faceLabels(mesh_.nFaces()/100);
 
     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
     {
-        label surfI = namedSurfaceIndex[faceI];
+        if (faceToZone[faceI] != -1)
-        if (surfI != -1)
         {
             // Free standing baffle?
             label ownZone = cellToZone[faceOwner[faceI]];
             label neiZone = cellToZone[faceNeighbour[faceI]];
-            if (max(ownZone, neiZone) == -1)
+            if (ownZone == neiZone)
             {
                 faceLabels.append(faceI);
             }
@@ -1843,13 +1808,12 @@ Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
         forAll(pp, i)
         {
             label faceI = pp.start()+i;
-            label surfI = namedSurfaceIndex[faceI];
+            if (faceToZone[faceI] != -1)
-            if (surfI != -1)
             {
                 // Free standing baffle?
                 label ownZone = cellToZone[faceOwner[faceI]];
                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
-                if (max(ownZone, neiZone) == -1)
+                if (ownZone == neiZone)
                 {
                     faceLabels.append(faceI);
                 }
@@ -2303,7 +2267,6 @@ void Foam::meshRefinement::consistentOrientation
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-// Split off unreachable areas of mesh.
 void Foam::meshRefinement::baffleAndSplitMesh
 (
     const bool doHandleSnapProblems,
@@ -2494,7 +2457,6 @@ void Foam::meshRefinement::baffleAndSplitMesh
 }
 
 
-// Split off (with optional buffer layers) unreachable areas of mesh.
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
 (
     const label nBufferLayers,
@@ -2778,8 +2740,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
 }
 
 
-// Find boundary points that connect to more than one cell region and
-// split them.
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints
 (
     const localPointRegion& regionSide
@@ -2833,8 +2793,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints
 }
 
 
-// Find boundary points that connect to more than one cell region and
-// split them.
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
 {
     // Analyse which points need to be duplicated
@@ -2844,7 +2802,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
 }
 
 
-// Zoning
 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 (
     const point& keepPoint,
@@ -3140,28 +3097,46 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
     }
 
+
+    //- Per face index of faceZone or -1
+    labelList faceToZone(mesh_.nFaces(), -1);
+
+    // Convert namedSurfaceIndex (index of named surfaces) to
+    // actual faceZone index
+
+    forAll(namedSurfaceIndex, faceI)
+    {
+        label surfI = namedSurfaceIndex[faceI];
+        if (surfI != -1)
+        {
+            faceToZone[faceI] = surfaceToFaceZone[surfI];
+        }
+    }
+
+
     // Topochange container
     polyTopoChange meshMod(mesh_);
 
 
 
     // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
-    labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
+    labelList neiCellZone;
+    syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
 
-        if (pp.coupled())
+        if (!pp.coupled())
         {
+            label bFaceI = pp.start()-mesh_.nInternalFaces();
             forAll(pp, i)
             {
-                label faceI = pp.start()+i;
+                neiCellZone[bFaceI++] = -1;
-                neiCellZone[faceI-mesh_.nInternalFaces()] =
-                    cellToZone[mesh_.faceOwner()[faceI]];
             }
         }
     }
-    syncTools::swapBoundaryFaceList(mesh_, neiCellZone);
+
+
 
     // Get per face whether is it master (of a coupled set of faces)
     const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
@@ -3188,7 +3163,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                 mesh_.faces(),
                 freeStandingBaffleFaces
                 (
-                    namedSurfaceIndex,
+                    faceToZone,
                     cellToZone,
                     neiCellZone
                 )
@@ -3202,18 +3177,25 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
             Info<< "Detected " << nFreeStanding << " free-standing zone faces"
                 << endl;
 
+            if (debug)
+            {
+                OBJstream str(mesh_.time().path()/"freeStanding.obj");
+                str.write(patch.localFaces(), patch.localPoints(), false);
+            }
+
+
             // Detect non-manifold edges
             labelList nMasterFacesPerEdge;
             calcPatchNumMasterFaces(isMasterFace, patch, nMasterFacesPerEdge);
 
             // Mark zones. Even a single original surface might create multiple
             // disconnected/non-manifold-connected zones
-            labelList faceToZone;
+            labelList faceToConnectedZone;
             const label nZones = markPatchZones
             (
                 patch,
                 nMasterFacesPerEdge,
-                faceToZone
+                faceToConnectedZone
             );
 
             Map<label> nPosOrientation(2*nZones);
@@ -3230,7 +3212,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                 isMasterFace,
                 patch,
                 nMasterFacesPerEdge,
-                faceToZone,
+                faceToConnectedZone,
                 nPosOrientation,
 
                 meshFlipMap
@@ -3254,7 +3236,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                         n = -1;
                     }
 
-                    nPosOrientation.find(faceToZone[faceI])() += n;
+                    nPosOrientation.find(faceToConnectedZone[faceI])() += n;
                 }
             }
             Pstream::mapCombineGather(nPosOrientation, plusEqOp<label>());
@@ -3281,7 +3263,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                 isMasterFace,
                 patch,
                 nMasterFacesPerEdge,
-                faceToZone,
+                faceToConnectedZone,
                 nPosOrientation,
 
                 meshFlipMap
@@ -3295,30 +3277,31 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
     for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
     {
-        label surfI = namedSurfaceIndex[faceI];
+        label faceZoneI = faceToZone[faceI];
 
-        if (surfI != -1)
+        if (faceZoneI != -1)
         {
             // Orient face zone to have slave cells in max cell zone.
+            // Note: logic to use flipMap should be consistent with logic
+            //       to pick up the freeStandingBaffleFaces!
+
             label ownZone = cellToZone[faceOwner[faceI]];
             label neiZone = cellToZone[faceNeighbour[faceI]];
 
             bool flip;
 
-            label maxZone = max(ownZone, neiZone);
+            if (ownZone == neiZone)
-
-            if (maxZone == -1)
             {
                 // free-standing face. Use geometrically derived orientation
                 flip = meshFlipMap[faceI];
             }
-            else if (ownZone == maxZone)
-            {
-                flip = false;
-            }
             else
             {
-                flip = true;
+                flip =
+                (
+                    ownZone == -1
+                 || (neiZone != -1 && ownZone > neiZone)
+                );
             }
 
             meshMod.setAction
@@ -3332,7 +3315,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                     false,                          // face flip
                     -1,                             // patch for face
                     false,                          // remove from zone
-                    surfaceToFaceZone[surfI],       // zone for face
+                    faceZoneI,                      // zone for face
                     flip                            // face flip in zone
                 )
             );
@@ -3349,35 +3332,27 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
         forAll(pp, i)
         {
-            label surfI = namedSurfaceIndex[faceI];
+            label faceZoneI = faceToZone[faceI];
 
-            if (surfI != -1)
+            if (faceZoneI != -1)
             {
                 label ownZone = cellToZone[faceOwner[faceI]];
                 label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
 
                 bool flip;
 
-                label maxZone = max(ownZone, neiZone);
+                if (ownZone == neiZone)
-
-                if (maxZone == -1)
                 {
                     // free-standing face. Use geometrically derived orientation
                     flip = meshFlipMap[faceI];
                 }
-                else if (ownZone == neiZone)
-                {
-                    // Free-standing zone face or coupled boundary. Keep master
-                    // face unflipped.
-                    flip = !isMasterFace[faceI];
-                }
-                else if (neiZone == maxZone)
-                {
-                    flip = true;
-                }
                 else
                 {
-                    flip = false;
+                    flip =
+                    (
+                        ownZone == -1
+                     || (neiZone != -1 && ownZone > neiZone)
+                    );
                 }
 
                 meshMod.setAction
@@ -3391,7 +3366,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                         false,                          // face flip
                         patchI,                         // patch for face
                         false,                          // remove from zone
-                        surfaceToFaceZone[surfI],       // zone for face
+                        faceZoneI,                      // zone for face
                         flip                            // face flip in zone
                     )
                 );