ENH: snappyHexMesh: allow cross-patch merging of patch faces. Fixes #1255.

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -1691,19 +1691,41 @@ int main(int argc, char *argv[])
     }
 
 
-    const bool mergePatchFaces
+    // How to treat co-planar faces
-    (
+    meshRefinement::FaceMergeType mergeType =
-        meshDict.lookupOrDefault("mergePatchFaces", true)
+        meshRefinement::FaceMergeType::GEOMETRIC;
-    );
-
-    if (!mergePatchFaces)
     {
-        Info<< "Not merging patch-faces of cell to preserve"
+        const bool mergePatchFaces
-            << " (split)hex cell shape."
+        (
-            << nl << endl;
+            meshDict.lookupOrDefault("mergePatchFaces", true)
+        );
+
+        if (!mergePatchFaces)
+        {
+            Info<< "Not merging patch-faces of cell to preserve"
+                << " (split)hex cell shape."
+                << nl << endl;
+            mergeType = meshRefinement::FaceMergeType::NONE;
+        }
+        else
+        {
+            const bool mergeAcrossPatches
+            (
+                meshDict.lookupOrDefault("mergeAcrossPatches", false)
+            );
+
+            if (mergeAcrossPatches)
+            {
+                Info<< "Merging co-planar patch-faces of cells"
+                    << ", regardless of patch assignment"
+                    << nl << endl;
+                mergeType = meshRefinement::FaceMergeType::IGNOREPATCH;
+            }
+        }
     }
 
 
+
     if (wantRefine)
     {
         cpuTime timer;
@@ -1732,7 +1754,7 @@ int main(int argc, char *argv[])
             refineParams,
             snapParams,
             refineParams.handleSnapProblems(),
-            mergePatchFaces,        // merge co-planar faces
+            mergeType,
             motionDict
         );
 
@@ -1784,7 +1806,7 @@ int main(int argc, char *argv[])
         (
             snapDict,
             motionDict,
-            mergePatchFaces,
+            mergeType,
             curvature,
             planarAngle,
             snapParams
@@ -1851,7 +1873,7 @@ int main(int argc, char *argv[])
             layerDict,
             motionDict,
             layerParams,
-            mergePatchFaces,
+            mergeType,
             preBalance,
             decomposer,
             distributor

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.C

@@ -129,6 +129,7 @@ bool Foam::combineFaces::validFace
 void Foam::combineFaces::regioniseFaces
 (
     const scalar minCos,
+    const bool mergeAcrossPatches,
     const label celli,
     const labelList& cEdges,
     Map<label>& faceRegion
@@ -143,16 +144,31 @@ void Foam::combineFaces::regioniseFaces
         label f0, f1;
         meshTools::getEdgeFaces(mesh_, celli, edgeI, f0, f1);
 
+        const vector& a0 = mesh_.faceAreas()[f0];
+        const vector& a1 = mesh_.faceAreas()[f1];
+
         const label p0 = patches.whichPatch(f0);
         const label p1 = patches.whichPatch(f1);
 
         // Face can be merged if
-        // - same non-coupled patch
         // - small angle
-        if (p0 != -1 && p0 == p1 && !patches[p0].coupled())
+        // - mergeAcrossPatches=false : same non-coupled patch
+        // - mergeAcrossPatches=true  : always
+        if
+        (
+            p0 != -1
+         && p1 != -1
+         && !patches[p0].coupled()
+         && !patches[p1].coupled()
+        )
         {
-            const vector f0Normal = normalised(mesh_.faceAreas()[f0]);
+            if (!mergeAcrossPatches && (p0 != p1))
-            const vector f1Normal = normalised(mesh_.faceAreas()[f1]);
+            {
+                continue;
+            }
+
+            const vector f0Normal = normalised(a0);
+            const vector f1Normal = normalised(a1);
 
             if ((f0Normal & f1Normal) > minCos)
             {
@@ -285,7 +301,8 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 (
     const scalar featureCos,
     const scalar minConcaveCos,
-    const labelHashSet& boundaryCells
+    const labelHashSet& boundaryCells,
+    const bool mergeAcrossPatches
 ) const
 {
     // Lists of faces that can be merged.
@@ -303,7 +320,14 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 
         // Region per face
         Map<label> faceRegion(cFaces.size());
-        regioniseFaces(featureCos, celli, cEdges, faceRegion);
+        regioniseFaces
+        (
+            featureCos,
+            mergeAcrossPatches,
+            celli,
+            cEdges,
+            faceRegion
+        );
 
         // Now we have in faceRegion for every face the region with planar
         // face sharing the same region. We now check whether the resulting
@@ -338,7 +362,7 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 
             // For every set check if it forms a valid convex face
 
-            forAllConstIters(regionToFaces, iter)
+            forAllIters(regionToFaces, iter)
             {
                 // Make face out of setFaces
                 indirectPrimitivePatch bigFace
@@ -354,7 +378,33 @@ Foam::labelListList Foam::combineFaces::getMergeSets
                 // Only store if -only one outside loop -which forms convex face
                 if (validFace(minConcaveCos, bigFace))
                 {
-                    allFaceSets.append(iter.val());
+                    labelList& faceIDs = iter.val();
+
+                    // For cross-patch merging we want to make the
+                    // largest face the one to decide the final patch
+                    // (i.e. master face)
+                    if (mergeAcrossPatches)
+                    {
+                        const vectorField& areas = mesh_.faceAreas();
+
+                        label maxIndex = 0;
+                        scalar maxMagSqr = magSqr(areas[faceIDs[0]]);
+                        for (label i = 1; i < faceIDs.size(); ++i)
+                        {
+                            const scalar a2 = magSqr(areas[faceIDs[i]]);
+                            if (a2 > maxMagSqr)
+                            {
+                                maxMagSqr = a2;
+                                maxIndex = i;
+                            }
+                        }
+                        if (maxIndex != 0)
+                        {
+                            Swap(faceIDs[0], faceIDs[maxIndex]);
+                        }
+                    }
+
+                    allFaceSets.append(faceIDs);
                 }
             }
         }
@@ -367,7 +417,8 @@ Foam::labelListList Foam::combineFaces::getMergeSets
 Foam::labelListList Foam::combineFaces::getMergeSets
 (
     const scalar featureCos,
-    const scalar minConcaveCos
+    const scalar minConcaveCos,
+    const bool mergeAcrossPatches
 ) const
 {
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
@@ -388,7 +439,13 @@ Foam::labelListList Foam::combineFaces::getMergeSets
         }
     }
 
-    return getMergeSets(featureCos, minConcaveCos, boundaryCells);
+    return getMergeSets
+    (
+        featureCos,
+        minConcaveCos,
+        boundaryCells,
+        mergeAcrossPatches
+    );
 }
 
 

src/dynamicMesh/polyTopoChange/polyTopoChange/combineFaces.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           |
+    \\  /    A nd           | Copyright (C)2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2016 OpenFOAM Foundation
@@ -103,6 +103,7 @@ class combineFaces
         void regioniseFaces
         (
             const scalar minCos,
+            const bool mergeAcrossPatches,
             const label celli,
             const labelList& cEdges,
             Map<label>& faceRegion
@@ -155,20 +156,27 @@ public:
         // Helper functions
 
             //- Extract lists of all (non-coupled) boundary faces on selected
-            //  cells that can be merged. Uses getFaceRegions.
+            //  cells that can be merged. Uses getFaceRegions. Optionally
+            //  allow faces-on-different-patches to be merged (into the largest
+            //  area face - could be improved). Note: causes a problem in
+            //  undoing - all restored faces get the patch/zone from the
+            //  master face.
             labelListList getMergeSets
             (
                 const scalar featureCos,
                 const scalar minConcaveCos,
-                const labelHashSet& boundaryCells
+                const labelHashSet& boundaryCells,
+                const bool mergeAcrossPatches = false
             ) const;
 
             //- Extract lists of all (non-coupled) boundary faces that can
-            //  be merged. Uses getFaceRegions.
+            //  be merged. Uses getFaceRegions. See note above about
+            //  mergeAcrossPatches.
             labelListList getMergeSets
             (
                 const scalar featureCos,
-                const scalar minConcaveCos
+                const scalar minConcaveCos,
+                const bool mergeAcrossPatches = false
             ) const;
 
             //- Gets outside of patch as a face (in mesh point labels)
@@ -197,7 +205,8 @@ public:
             //  Returns maps from added restored point to
             //  original point label (i.e. content of savedPointLabels_).
             //  (only restoredPoints are actually set; rest are just for
-            //   generalness)
+            //   generalness). See note above about restoring faces from
+            //  different patches (mergeAcrossPatches)
             void setUnrefinement
             (
                 const labelList& masterFaces,

src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2015-2019 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2015-2017 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2017 OpenFOAM Foundation
@@ -129,6 +129,16 @@ public:
             REMOVE = 4      //!< set value to -1 any face that was refined
         };
 
+        //- Enumeration for what to do with co-planar patch faces on a single
+        //  cell
+        enum FaceMergeType
+        {
+            NONE,           // no merging
+            GEOMETRIC,      // use feature angle
+            IGNOREPATCH     // use feature angle, allow merging of different
+                            // patches
+        };
+
 
 private:
 
@@ -1463,7 +1473,8 @@ public:
                     const scalar minCos,
                     const scalar concaveCos,
                     const label mergeSize,
-                    const labelList& patchIDs
+                    const labelList& patchIDs,
+                    const meshRefinement::FaceMergeType mergeType
                 );
 
                 //- Merge coplanar faces. preserveFaces is != -1 for faces
@@ -1474,7 +1485,8 @@ public:
                     const scalar concaveCos,
                     const labelList& patchIDs,
                     const dictionary& motionDict,
-                    const labelList& preserveFaces
+                    const labelList& preserveFaces,
+                    const meshRefinement::FaceMergeType mergeType
                 );
 
                 autoPtr<mapPolyMesh> doRemovePoints

src/mesh/snappyHexMesh/meshRefinement/meshRefinementMerge.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2016-2018 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2016-2019 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
                             | Copyright (C) 2011-2014 OpenFOAM Foundation
@@ -42,7 +42,8 @@ Foam::label Foam::meshRefinement::mergePatchFaces
     const scalar minCos,
     const scalar concaveCos,
     const label mergeSize,
-    const labelList& patchIDs
+    const labelList& patchIDs,
+    const meshRefinement::FaceMergeType mergeType
 )
 {
     // Patch face merging engine
@@ -73,7 +74,8 @@ Foam::label Foam::meshRefinement::mergePatchFaces
         (
             minCos,
             concaveCos,
-            boundaryCells
+            boundaryCells,
+            (mergeType == FaceMergeType::IGNOREPATCH) // merge across patches?
         )
     );
 
@@ -249,7 +251,8 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
     const scalar concaveCos,
     const labelList& patchIDs,
     const dictionary& motionDict,
-    const labelList& preserveFaces
+    const labelList& preserveFaces,
+    const meshRefinement::FaceMergeType mergeType
 )
 {
     // Patch face merging engine
@@ -286,7 +289,8 @@ Foam::label Foam::meshRefinement::mergePatchFacesUndo
         (
             minCos,
             concaveCos,
-            boundaryCells
+            boundaryCells,
+            (mergeType == FaceMergeType::IGNOREPATCH) // merge across patches?
         )
     );
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.C

@@ -3219,7 +3219,8 @@ Foam::snappyLayerDriver::snappyLayerDriver
 void Foam::snappyLayerDriver::mergePatchFacesUndo
 (
     const layerParameters& layerParams,
-    const dictionary& motionDict
+    const dictionary& motionDict,
+    const meshRefinement::FaceMergeType mergeType
 )
 {
     // Clip to 30 degrees. Not helpful!
@@ -3260,7 +3261,8 @@ void Foam::snappyLayerDriver::mergePatchFacesUndo
         concaveCos,
         meshRefiner_.meshedPatches(),
         motionDict,
-        duplicateFace
+        duplicateFace,
+        mergeType   // How to merge co-planar patch faces
     );
 
     nChanged += meshRefiner_.mergeEdgesUndo(minCos, motionDict);
@@ -4635,7 +4637,7 @@ void Foam::snappyLayerDriver::doLayers
     const dictionary& shrinkDict,
     const dictionary& motionDict,
     const layerParameters& layerParams,
-    const bool mergePatchFaces,
+    const meshRefinement::FaceMergeType mergeType,
     const bool preBalance,
     decompositionMethod& decomposer,
     fvMeshDistribute& distributor
@@ -4653,9 +4655,13 @@ void Foam::snappyLayerDriver::doLayers
     Info<< "Using mesh parameters " << motionDict << nl << endl;
 
     // Merge coplanar boundary faces
-    if (mergePatchFaces)
+    if
+    (
+        mergeType == meshRefinement::FaceMergeType::GEOMETRIC
+     || mergeType == meshRefinement::FaceMergeType::IGNOREPATCH
+    )
     {
-        mergePatchFacesUndo(layerParams, motionDict);
+        mergePatchFacesUndo(layerParams, motionDict, mergeType);
     }
 
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyLayerDriver.H

@@ -634,7 +634,8 @@ public:
             void mergePatchFacesUndo
             (
                 const layerParameters& layerParams,
-                const dictionary& motionDict
+                const dictionary& motionDict,
+                const meshRefinement::FaceMergeType mergeType
             );
 
             //- Add cell layers
@@ -654,12 +655,11 @@ public:
                 const dictionary& shrinkDict,
                 const dictionary& motionDict,
                 const layerParameters& layerParams,
-                const bool mergePatchFaces,         // merging patch faces
+                const meshRefinement::FaceMergeType mergeType,
                 const bool preBalance,              // balance before adding?
                 decompositionMethod& decomposer,
                 fvMeshDistribute& distributor
             );
-
 };
 
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C

@@ -2794,7 +2794,7 @@ void Foam::snappyRefineDriver::addFaceZones
 
 void Foam::snappyRefineDriver::mergePatchFaces
 (
-    const bool geometricMerge,
+    const meshRefinement::FaceMergeType mergeType,
     const refinementParameters& refineParams,
     const dictionary& motionDict
 )
@@ -2812,7 +2812,11 @@ void Foam::snappyRefineDriver::mergePatchFaces
 
     const fvMesh& mesh = meshRefiner_.mesh();
 
-    if (geometricMerge)
+    if
+    (
+        mergeType == meshRefinement::FaceMergeType::GEOMETRIC
+     || mergeType == meshRefinement::FaceMergeType::IGNOREPATCH
+    )
     {
         meshRefiner_.mergePatchFacesUndo
         (
@@ -2820,7 +2824,8 @@ void Foam::snappyRefineDriver::mergePatchFaces
             Foam::cos(degToRad(45.0)),
             meshRefiner_.meshedPatches(),
             motionDict,
-            labelList(mesh.nFaces(), -1)
+            labelList(mesh.nFaces(), -1),
+            mergeType
         );
     }
     else
@@ -2831,7 +2836,8 @@ void Foam::snappyRefineDriver::mergePatchFaces
             Foam::cos(degToRad(45.0)),
             Foam::cos(degToRad(45.0)),
             4,          // only merge faces split into 4
-            meshRefiner_.meshedPatches()
+            meshRefiner_.meshedPatches(),
+            meshRefinement::FaceMergeType::GEOMETRIC // no merge across patches
         );
     }
 
@@ -2855,7 +2861,7 @@ void Foam::snappyRefineDriver::doRefine
     const refinementParameters& refineParams,
     const snapParameters& snapParams,
     const bool prepareForSnapping,
-    const bool doMergePatchFaces,
+    const meshRefinement::FaceMergeType mergeType,
     const dictionary& motionDict
 )
 {
@@ -3063,7 +3069,7 @@ void Foam::snappyRefineDriver::doRefine
     // Do something about cells with refined faces on the boundary
     if (prepareForSnapping)
     {
-        mergePatchFaces(doMergePatchFaces, refineParams, motionDict);
+        mergePatchFaces(mergeType, refineParams, motionDict);
     }
 
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.H

@@ -43,6 +43,7 @@ SourceFiles
 #include "writer.H"
 #include "DynamicList.H"
 #include "labelVector.H"
+#include "meshRefinement.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -53,7 +54,6 @@ namespace Foam
 class refinementParameters;
 class snapParameters;
 
-class meshRefinement;
 class decompositionMethod;
 class fvMeshDistribute;
 class fvMesh;
@@ -226,7 +226,7 @@ class snappyRefineDriver
         //- Merge refined boundary faces (from exposing coarser cell)
         void mergePatchFaces
         (
-            const bool geometricMerge,
+            const meshRefinement::FaceMergeType mergeType,
             const refinementParameters& refineParams,
             const dictionary& motionDict
         );
@@ -268,7 +268,7 @@ public:
             const refinementParameters& refineParams,
             const snapParameters& snapParams,
             const bool prepareForSnapping,
-            const bool mergePatchFaces,
+            const meshRefinement::FaceMergeType mergeType,
             const dictionary& motionDict
         );
 

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.C

@@ -2528,7 +2528,7 @@ void Foam::snappySnapDriver::doSnap
 (
     const dictionary& snapDict,
     const dictionary& motionDict,
-    const bool mergePatchFaces,
+    const meshRefinement::FaceMergeType mergeType,
     const scalar featureCos,
     const scalar planarAngle,
     const snapParameters& snapParams
@@ -3031,7 +3031,11 @@ void Foam::snappySnapDriver::doSnap
         repatchToSurface(snapParams, adaptPatchIDs, duplicateFace);
     }
 
-    if (mergePatchFaces)
+    if
+    (
+        mergeType == meshRefinement::FaceMergeType::GEOMETRIC
+     || mergeType == meshRefinement::FaceMergeType::IGNOREPATCH
+    )
     {
         labelList duplicateFace(getInternalOrBaffleDuplicateFace());
 
@@ -3044,7 +3048,8 @@ void Foam::snappySnapDriver::doSnap
             featureCos,     // concaveCos
             meshRefiner_.meshedPatches(),
             motionDict,
-            duplicateFace   // faces not to merge
+            duplicateFace,  // faces not to merge
+            mergeType
         );
 
         nChanged += meshRefiner_.mergeEdgesUndo(featureCos, motionDict);

src/mesh/snappyHexMesh/snappyHexMeshDriver/snappySnapDriver.H

@@ -781,12 +781,11 @@ public:
             (
                 const dictionary& snapDict,
                 const dictionary& motionDict,
-                const bool mergePatchFaces,
+                const meshRefinement::FaceMergeType mergeType,
                 const scalar featureCos,
                 const scalar planarAngle,
                 const snapParameters& snapParams
             );
-
 };