From bf43a9d91eb3fd8ac05ca1b268a4dc2111d1977b Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 30 Aug 2018 11:38:38 +0100 Subject: [PATCH] ENH: snappyHexMesh: allow medialAxis to shrink 1 or 2 cell thick channel. Fixes #986. --- .../snappyHexMesh/snappyHexMeshDict | 7 ++ .../medialAxisMeshMover.C | 114 ++++++++++++------ .../medialAxisMeshMover.H | 3 +- 3 files changed, 86 insertions(+), 38 deletions(-) diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index 8d521b254c..cef09c6e94 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -665,6 +665,13 @@ addLayersControls // Default is 0.5*featureAngle. Set to -180 always attempt extrusion //layerTerminationAngle 25; + // Optional: disable shrinking of edges that have one (or two) points + // on an extruded patch. + // Default is false to enable single/two cell thick channels to still + // have layers. In <=1806 this was true by default. On larger gaps it + // should have no effect. + //disableWallEdges true; + // Optional: at non-patched sides allow mesh to slip if extrusion // direction makes angle larger than slipFeatureAngle. Default is // 0.5*featureAngle. diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C index 40095d5233..6bb99d604d 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.C @@ -58,33 +58,42 @@ bool Foam::medialAxisMeshMover::isMaxEdge ( const List>& pointWallDist, const label edgeI, - const scalar minCos + const scalar minCos, + const bool disableWallEdges ) const { const pointField& points = mesh().points(); - - // Do not mark edges with one side on moving wall. - const edge& e = mesh().edges()[edgeI]; - vector v0(points[e[0]] - pointWallDist[e[0]].origin()); - scalar magV0(mag(v0)); - - if (magV0 < SMALL) + if (disableWallEdges) { - return false; + // 1. Do not mark edges with one side on moving wall. + vector v0(points[e[0]] - pointWallDist[e[0]].origin()); + scalar magV0(mag(v0)); + if (magV0 < SMALL) + { + return false; + } + + vector v1(points[e[1]] - pointWallDist[e[1]].origin()); + scalar magV1(mag(v1)); + if (magV1 < SMALL) + { + return false; + } } - vector v1(points[e[1]] - pointWallDist[e[1]].origin()); - scalar magV1(mag(v1)); + //// 2. Do not mark edges with both sides on a moving wall. + //vector v0(points[e[0]] - pointWallDist[e[0]].origin()); + //scalar magV0(mag(v0)); + //vector v1(points[e[1]] - pointWallDist[e[1]].origin()); + //scalar magV1(mag(v1)); + //if (magV0 < SMALL && magV1 < SMALL) + //{ + // return false; + //} - if (magV1 < SMALL) - { - return false; - } - - - //- Detect based on vector to nearest point differing for both endpoints + //// 3. Detect based on vector to nearest point differing for both endpoints //v0 /= magV0; //v1 /= magV1; // @@ -170,6 +179,13 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict) mesh().globalData().nTotalPoints() ); + const bool disableWallEdges = coeffDict.lookupOrDefault + ( + "disableWallEdges", + false + ); + + // Predetermine mesh edges // ~~~~~~~~~~~~~~~~~~~~~~~ @@ -336,7 +352,16 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict) } } - else if (isMaxEdge(pointWallDist, edgeI, minMedialAxisAngleCos)) + else if + ( + isMaxEdge + ( + pointWallDist, + edgeI, + minMedialAxisAngleCos, + disableWallEdges + ) + ) { // Both end points of edge have very different nearest wall // point. Mark both points as medial axis points. @@ -351,41 +376,56 @@ void Foam::medialAxisMeshMover::update(const dictionary& coeffDict) // Calculate distance along edge const point& p0 = points[e[0]]; + const point& origin0 = pointWallDist[e[0]].origin(); const point& p1 = points[e[1]]; - scalar dist0 = (p0-pointWallDist[e[0]].origin()) & eVec; - scalar dist1 = (pointWallDist[e[1]].origin()-p1) & eVec; + const point& origin1 = pointWallDist[e[1]].origin(); + scalar dist0 = (p0-origin0) & eVec; + scalar dist1 = (origin1-p1) & eVec; scalar s = 0.5*(dist1+eMag+dist0); - point medialAxisPt; + point medialAxisPt(vector::max); if (s <= dist0) { - medialAxisPt = p0; + // Make sure point is not on wall. Note that this + // check used to be inside isMaxEdge. + if (magSqr((p0-origin0)) > Foam::sqr(SMALL)) + { + medialAxisPt = p0; + } } else if (s >= dist0+eMag) { - medialAxisPt = p1; + // Make sure point is not on wall. Note that this + // check used to be inside isMaxEdge. + if (magSqr((p1-origin1)) > Foam::sqr(SMALL)) + { + medialAxisPt = p1; + } } else { medialAxisPt = p0+(s-dist0)*eVec; } - forAll(e, ep) + if (medialAxisPt != vector::max) { - label pointI = e[ep]; - - if (!pointMedialDist[pointI].valid(dummyTrackData)) + forAll(e, ep) { - maxPoints.append(pointI); - maxInfo.append - ( - pointEdgePoint + label pointI = e[ep]; + + if (!pointMedialDist[pointI].valid(dummyTrackData)) + { + maxPoints.append(pointI); + maxInfo.append ( - medialAxisPt, //points[pointI], - magSqr(points[pointI]-medialAxisPt)//0.0 - ) - ); - pointMedialDist[pointI] = maxInfo.last(); + pointEdgePoint + ( + medialAxisPt, //points[pointI], + magSqr(points[pointI]-medialAxisPt)//0.0 + ) + ); + pointMedialDist[pointI] = maxInfo.last(); + } } } } diff --git a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H index 6fa9a02405..8f3ea990f4 100644 --- a/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H +++ b/src/mesh/snappyHexMesh/externalDisplacementMeshMover/medialAxisMeshMover.H @@ -114,7 +114,8 @@ class medialAxisMeshMover ( const List>& pointWallDist, const label edgeI, - const scalar minCos + const scalar minCos, + const bool checkWall ) const; //- Initialise medial axis. Uses pointDisplacement field only