From 25a373bd395d1e8335849cd159a1f29c7f7d216f Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 16 May 2013 12:49:16 +0100 Subject: [PATCH] ENH: snappyHexMesh: enable free-standing baffle merging --- .../snappyHexMesh/snappyHexMeshDict | 6 + .../autoHexMeshDriver/autoRefineDriver.C | 5 +- .../refinementParameters.C | 11 +- .../refinementParameters.H | 11 +- .../meshRefinement/meshRefinement.H | 35 +- .../meshRefinement/meshRefinementBaffles.C | 612 +++++++++++------- .../meshRefinement/meshRefinementRefine.C | 25 +- 7 files changed, 445 insertions(+), 260 deletions(-) diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index f448faa40c..ca4a9dc6b5 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -180,6 +180,12 @@ castellatedMeshControls // - used if feature snapping (see snapControls below) is used resolveFeatureAngle 30; + // Planar angle: + // - used to determine if neighbouring surface normals + // are roughly the same so e.g. free-standing baffles can be merged + // If not specified same as resolveFeatureAngle + //planarAngle 30; + // Region-wise refinement // ~~~~~~~~~~~~~~~~~~~~~~ diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C index 484f851c0f..75ea461102 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C @@ -522,6 +522,7 @@ void Foam::autoRefineDriver::baffleAndSplitMesh false, // perpendicular edge connected cells scalarField(0), // per region perpendicular angle !handleSnapProblems, // merge free standing baffles? + refineParams.planarAngle(), motionDict, const_cast(mesh.time()), globalToMasterPatch_, @@ -606,8 +607,8 @@ void Foam::autoRefineDriver::splitAndMergeBaffles handleSnapProblems, handleSnapProblems, // remove perp edge connected cells perpAngle, // perp angle - false, // merge free standing baffles? - //true, // merge free standing baffles? + true, // merge free standing baffles? + refineParams.planarAngle(), // planar angle motionDict, const_cast(mesh.time()), globalToMasterPatch_, diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C index ab2140a4ec..7bf8d8cf6c 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -41,6 +41,7 @@ Foam::refinementParameters::refinementParameters maxLocalCells_(readLabel(dict.lookup("procCellLimit"))), minRefineCells_(readLabel(dict.lookup("minimumRefine"))), curvature_(readScalar(dict.lookup("curvature"))), + planarAngle_(dict.lookupOrDefault("planarAngle", curvature_)), nBufferLayers_(readLabel(dict.lookup("nBufferLayers"))), keepPoints_(dict.lookup("keepPoints")), allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")), @@ -53,6 +54,14 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict) maxGlobalCells_(readLabel(dict.lookup("maxGlobalCells"))), maxLocalCells_(readLabel(dict.lookup("maxLocalCells"))), minRefineCells_(readLabel(dict.lookup("minRefinementCells"))), + planarAngle_ + ( + dict.lookupOrDefault + ( + "planarAngle", + readScalar(dict.lookup("resolveFeatureAngle")) + ) + ), nBufferLayers_(readLabel(dict.lookup("nCellsBetweenLevels"))), keepPoints_(pointField(1, dict.lookup("locationInMesh"))), allowFreeStandingZoneFaces_(dict.lookup("allowFreeStandingZoneFaces")), diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H index 66867e7f55..4fb05644e8 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -67,6 +67,9 @@ class refinementParameters //- Curvature scalar curvature_; + //- Planarity criterion + scalar planarAngle_; + //- Number of layers between different refinement levels const label nBufferLayers_; @@ -129,6 +132,12 @@ public: return curvature_; } + //- Angle when two intersections are considered to be planar + scalar planarAngle() const + { + return planarAngle_; + } + //- Number of layers between different refinement levels label nBufferLayers() const { diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H index 6ab861d8a4..60e6261ccc 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H @@ -253,6 +253,14 @@ private: label& nRefine ); + //- Mark every cell with level of feature passing through it + // (or -1 if not passed through). Uses tracking. + void markFeatureCellLevel + ( + const point& keepPoint, + labelList& maxFeatureLevel + ) const; + //- Calculate list of cells to refine based on intersection of // features. label markFeatureRefinement @@ -337,16 +345,6 @@ private: const labelList& globalToSlavePatch ) const; - //- Geometric test on see whether face needs to be baffled: - // is face boundary face and perpendicular to surface normal? - bool validBaffleTopology - ( - const label faceI, - const vector& n1, - const vector& n2, - const vector& testDir - ) const; - //- Determine patches for baffles void getBafflePatches ( @@ -435,7 +433,11 @@ private: //- Extract those baffles (duplicate) faces that are on the edge // of a baffle region. These are candidates for merging. - List freeStandingBaffles(const List&) const; + List freeStandingBaffles + ( + const List&, + const scalar freeStandingAngle + ) const; // Zone handling @@ -489,6 +491,16 @@ private: labelList& namedSurfaceIndex ) const; + //- Remove any loose standing cells + void handleSnapProblems + ( + const bool removeEdgeConnectedCells, + const scalarField& perpendicularAngle, + const dictionary& motionDict, + Time& runTime, + const labelList& globalToMasterPatch, + const labelList& globalToSlavePatch + ); //- Disallow default bitwise copy construct meshRefinement(const meshRefinement&); @@ -707,6 +719,7 @@ public: const bool removeEdgeConnectedCells, const scalarField& perpendicularAngle, const bool mergeFreeStanding, + const scalar freeStandingAngle, const dictionary& motionDict, Time& runTime, const labelList& globalToMasterPatch, diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C index 339e51cd7f..7a4f87753d 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C @@ -44,6 +44,7 @@ License #include "regionSplit.H" #include "removeCells.H" #include "unitConversion.H" +#include "OBJstream.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -217,43 +218,43 @@ Foam::label Foam::meshRefinement::getBafflePatch } -// 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; -} +//// 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 @@ -800,7 +801,8 @@ Foam::autoPtr Foam::meshRefinement::createZoneBaffles // region. Foam::List Foam::meshRefinement::freeStandingBaffles ( - const List& couples + const List& couples, + const scalar planarAngle ) const { // All duplicate faces on edge of the patch are to be merged. @@ -809,6 +811,22 @@ Foam::List Foam::meshRefinement::freeStandingBaffles labelList nBafflesPerEdge(mesh_.nEdges(), 0); + // This algorithm is quite tricky. We don't want to use edgeFaces and + // also want it to run in parallel so it is now an algorithm over + // all (boundary) faces instead. + // We want to pick up any edges that are only used by the baffle + // or internal faces but not by any other boundary faces. So + // - increment count on an edge by 1 if it is used by any (uncoupled) + // boundary face. + // - increment count on an edge by 1000000 if it is used by a baffle face + // - sum in parallel + // + // So now any edge that is used by baffle faces only will have the + // value 2*1000000+2*1. + + + const label baffleValue = 1000000; + // Count number of boundary faces per edge // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -847,18 +865,22 @@ Foam::List Foam::meshRefinement::freeStandingBaffles forAll(couples, i) { - const labelList& fEdges0 = mesh_.faceEdges(couples[i].first(), fe0); - - forAll(fEdges0, fEdgeI) { - nBafflesPerEdge[fEdges0[fEdgeI]] += 1000000; + label f0 = couples[i].first(); + const labelList& fEdges0 = mesh_.faceEdges(f0, fe0); + forAll(fEdges0, fEdgeI) + { + nBafflesPerEdge[fEdges0[fEdgeI]] += baffleValue; + } } - const labelList& fEdges1 = mesh_.faceEdges(couples[i].second(), fe1); - - forAll(fEdges1, fEdgeI) { - nBafflesPerEdge[fEdges1[fEdgeI]] += 1000000; + label f1 = couples[i].second(); + const labelList& fEdges1 = mesh_.faceEdges(f1, fe1); + forAll(fEdges1, fEdgeI) + { + nBafflesPerEdge[fEdges1[fEdgeI]] += baffleValue; + } } } @@ -873,8 +895,8 @@ Foam::List Foam::meshRefinement::freeStandingBaffles // Baffles which are not next to other boundaries and baffles will have - // nBafflesPerEdge value 2*1000000+2*1 (from 2 boundary faces which are - // both baffle faces) + // nBafflesPerEdge value 2*baffleValue+2*1 (from 2 boundary faces which + // are both baffle faces) List filteredCouples(couples.size()); label filterI = 0; @@ -889,15 +911,15 @@ Foam::List Foam::meshRefinement::freeStandingBaffles == patches.whichPatch(couple.second()) ) { - const labelList& fEdges = mesh_.faceEdges(couples[i].first()); + const labelList& fEdges = mesh_.faceEdges(couple.first()); forAll(fEdges, fEdgeI) { label edgeI = fEdges[fEdgeI]; - if (nBafflesPerEdge[edgeI] == 2*1000000+2*1) + if (nBafflesPerEdge[edgeI] == 2*baffleValue+2*1) { - filteredCouples[filterI++] = couples[i]; + filteredCouples[filterI++] = couple; break; } } @@ -905,111 +927,127 @@ Foam::List Foam::meshRefinement::freeStandingBaffles } filteredCouples.setSize(filterI); - //Info<< "freeStandingBaffles : from " - // << returnReduce(couples.size(), sumOp