From 16b7707c31200ac0b6276cd0a58ffaa3a39c75ad Mon Sep 17 00:00:00 2001 From: mattijs Date: Thu, 4 Oct 2018 12:03:53 +0100 Subject: [PATCH] ENH: snappyHexMesh: directional smoothing. See #1031 --- .../shellSurfaces/shellSurfaces.C | 19 +- .../shellSurfaces/shellSurfaces.H | 6 + .../snappyHexMeshDriver/snappyRefineDriver.C | 832 ++++++++---------- .../snappyHexMeshDriver/snappyRefineDriver.H | 25 +- .../Allrun | 22 +- .../system/decomposeParDict | 9 - .../system/snappyHexMeshDict | 16 +- 7 files changed, 434 insertions(+), 495 deletions(-) diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C index 81e52c58f6..ff0f06353c 100644 --- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C +++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2015-2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2015-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -587,6 +587,7 @@ Foam::shellSurfaces::shellSurfaces dirLevels_.setSize(shellI); smoothDirection_.setSize(shellI); nSmoothExpansion_.setSize(shellI); + nSmoothPosition_.setSize(shellI); extendedGapLevel_.setSize(shellI); extendedGapMode_.setSize(shellI); @@ -666,10 +667,16 @@ Foam::shellSurfaces::shellSurfaces // Directional smoothing // ~~~~~~~~~~~~~~~~~~~~~ + nSmoothExpansion_[shellI] = 0; + nSmoothPosition_[shellI] = 0; smoothDirection_[shellI] = dict.lookupOrDefault("smoothDirection", vector::zero); - nSmoothExpansion_[shellI] = - dict.lookupOrDefault("nSmoothExpansion", 0); + + if (smoothDirection_[shellI] != vector::zero) + { + dict.lookup("nSmoothExpansion") >> nSmoothExpansion_[shellI]; + dict.lookup("nSmoothPosition") >> nSmoothPosition_[shellI]; + } // Gap specification @@ -824,6 +831,12 @@ const Foam::vectorField& Foam::shellSurfaces::smoothDirection() const } +const Foam::labelList& Foam::shellSurfaces::nSmoothPosition() const +{ + return nSmoothPosition_; +} + + void Foam::shellSurfaces::findHigherLevel ( const pointField& pt, diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H index 8b4f8a8ef9..2390a98e95 100644 --- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H +++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.H @@ -98,6 +98,9 @@ private: //- Per shell the directional smoothing iterations labelList nSmoothExpansion_; + //- Per shell the positional smoothing iterations + labelList nSmoothPosition_; + // Gap level refinement @@ -246,6 +249,9 @@ public: //- Per shell the directional smoothing iterations const labelList& nSmoothExpansion() const; + + //- Per shell the positional smoothing iterations + const labelList& nSmoothPosition() const; }; diff --git a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C index 21ca736aa3..c163f39eda 100644 --- a/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C +++ b/src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C @@ -41,17 +41,14 @@ License #include "labelVector.H" #include "profiling.H" #include "searchableSurfaces.H" -#include "PointIntegrateData.H" -#include "PointEdgeWave.H" -#include "OBJstream.H" +#include "fvMeshSubset.H" +#include "interpolationTable.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - -defineTypeNameAndDebug(snappyRefineDriver, 0); - + defineTypeNameAndDebug(snappyRefineDriver, 0); } // End namespace Foam @@ -1788,81 +1785,102 @@ Foam::label Foam::snappyRefineDriver::directionalShellRefine } -void Foam::snappyRefineDriver::getVirtualEdgeLength +void Foam::snappyRefineDriver::mergeAndSmoothRatio ( - const vector& userDirection, - const labelList& pointLabels, - const PackedBoolList& isXFace, - scalarField& maxUserSize -) const + const scalarList& allSeedPointDist, + const label nSmoothExpansion, + List>& keyAndValue +) { - const fvMesh& mesh = meshRefiner_.mesh(); - const labelList& own = mesh.faceOwner(); - const labelList& nei = mesh.faceNeighbour(); - const vectorField& fA = mesh.faceAreas(); - const pointField& faceCentres = mesh.faceCentres(); - const pointField& cellCentres = mesh.cellCentres(); + // Merge duplicate distance from coupled locations to get unique + // distances to operate on, do on master + SortableList unmergedDist(allSeedPointDist); + DynamicList mergedDist; - scalarField sumWeights(mesh.nCells(), 0.0); - scalarField sumArea(mesh.nCells(), 0.0); - - forAll(own, facei) + scalar prevDist = GREAT; + forAll(unmergedDist, i) { - const scalar fc = (faceCentres[facei]&userDirection); - const scalar cc = (cellCentres[own[facei]]&userDirection); - // Calculate 'pyramid volume' - scalar pyr3Vol = (fA[facei]&userDirection)*(fc-cc); - sumWeights[own[facei]] += pyr3Vol; - sumArea[own[facei]] += mag(fA[facei]&userDirection); - } - - forAll(nei, facei) - { - const scalar fc = (faceCentres[facei]&userDirection); - const scalar cc = (cellCentres[nei[facei]]&userDirection); - // Calculate 'pyramid volume' - scalar pyr3Vol = (fA[facei]&userDirection)*(cc-fc); - sumWeights[nei[facei]] += pyr3Vol; - sumArea[nei[facei]] += mag(fA[facei]&userDirection); - } - - const scalarField cellSize(2*sumWeights/sumArea); - - maxUserSize.setSize(mesh.nPoints()); - maxUserSize = -GREAT; - - forAll(pointLabels, i) - { - label pointi = pointLabels[i]; - const labelList& pFaces = mesh.pointFaces()[pointi]; - scalar& maxSz = maxUserSize[pointi]; - forAll(pFaces, pFacei) + scalar curDist = unmergedDist[i]; + scalar difference = mag(curDist - prevDist); + if (difference > meshRefiner_.mergeDistance()) + //if (difference > 0.01) { - label facei = pFaces[pFacei]; - if (isXFace[facei]) - { - const point& fc = faceCentres[facei]; - if (((cellCentres[own[facei]]-fc)&userDirection) > 0.0) - { - maxSz = max(maxSz, cellSize[own[facei]]); - } - if (mesh.isInternalFace(facei)) - { - if (((cellCentres[nei[facei]]-fc)&userDirection) > 0.0) - { - maxSz = max(maxSz, cellSize[nei[facei]]); - } - } - } + mergedDist.append(curDist); + prevDist = curDist; } } - syncTools::syncPointList - ( - mesh, - maxUserSize, - maxEqOp(), - -GREAT - ); + + // Sort the unique distances + SortableList sortedDist(mergedDist); + labelList indexSet = sortedDist.indices(); + + // Get updated position starting from original (undistorted) mesh + scalarList seedPointsNewLocation = sortedDist; + + scalar initResidual = 0.0; + scalar prevIterResidual = GREAT; + + for (label iter = 0; iter < nSmoothExpansion; iter++) + { + + // Position based edge averaging algorithm operated on + // all seed plane locations in normalized form. + // + // 0 1 2 3 4 5 6 (edge numbers) + // ---x---x---x---x---x---x--- + // 0 1 2 3 4 5 (point numbers) + // + // Average of edge 1-3 in terms of position + // = (point3 - point0)/3 + // Keeping points 0-1 frozen, new position of point 2 + // = position2 + (average of edge 1-3 as above) + for(label i = 2; i 0) + if + ( + shells.nSmoothExpansion()[shellI] > 0 + || shells.nSmoothPosition()[shellI] > 0 + ) { label surfi = shells.shells()[shellI]; const vector& userDirection = shells.smoothDirection()[shellI]; @@ -1894,33 +1915,48 @@ Foam::label Foam::snappyRefineDriver::directionalSmooth // Extract inside points labelList pointLabels; - PackedBoolList isInsidePoint(mesh.nPoints()); { // Get inside points List volType; - geometry[surfi].getVolumeType(points, volType); + geometry[surfi].getVolumeType(baseMesh.points(), volType); + label nInside = 0; forAll(volType, pointi) { if (volType[pointi] == volumeType::INSIDE) { - isInsidePoint.set(pointi); + nInside++; } } - pointLabels = isInsidePoint.used(); + pointLabels.setSize(nInside); + nInside = 0; + forAll(volType, pointi) + { + if (volType[pointi] == volumeType::INSIDE) + { + pointLabels[nInside++] = pointi; + } + } + + //PackedBoolList isInsidePoint(baseMesh.nPoints()); + //forAll(volType, pointi) + //{ + // if (volType[pointi] == volumeType::INSIDE) + // { + // isInsidePoint.set(pointi); + // } + //} + //pointLabels = isInsidePoint.used(); } - - const edgeList& edges = mesh.edges(); - const labelListList& pointEdges = mesh.pointEdges(); - // Mark all directed edges - PackedBoolList isXEdge(edges.size()); + PackedBoolList isXEdge(baseMesh.edges().size()); { + const edgeList& edges = baseMesh.edges(); forAll(edges, edgei) { const edge& e = edges[edgei]; - vector eVec(e.vec(points)); + vector eVec(e.vec(baseMesh.points())); eVec /= mag(eVec); if (mag(eVec&userDirection) > 0.9) { @@ -1929,343 +1965,190 @@ Foam::label Foam::snappyRefineDriver::directionalSmooth } } + // Get the extreme of smoothing region and + // normalize all points within + const scalar totalLength = + geometry[surfi].bounds().span() + & userDirection; + const scalar startPosition = + geometry[surfi].bounds().min() + & userDirection; - // Mark all directed faces (with all points inside the region) - PackedBoolList isXFace(mesh.nFaces()); - { - const vectorField& faceAreas = mesh.faceAreas(); - const faceList& faces = mesh.faces(); - - forAll(faces, facei) - { - const face& f = faces[facei]; - - bool allInRegion = true; - for (label pointi : f) - { - if (!isInsidePoint[pointi]) - { - allInRegion = false; - break; - } - } - - if (allInRegion) - { - vector n(faceAreas[facei]); - n /= mag(n); - if (mag(n&userDirection) > 0.9) - { - isXFace[facei] = true; - } - } - } - } - - - // Get left and right local edge - labelList leftEdge(points.size(), -1); - labelList rightEdge(points.size(), -1); - + scalarField normalizedPosition(pointLabels.size(), 0); forAll(pointLabels, i) { label pointi = pointLabels[i]; - const labelList& pEdges = pointEdges[pointi]; - forAll(pEdges, pEdgei) - { - label edgei = pEdges[pEdgei]; - const edge& e = edges[edgei]; - label otherPointi = e.otherVertex(pointi); - vector eVec(points[otherPointi]-points[pointi]); - eVec /= mag(eVec); - scalar cosAngle = (eVec&userDirection); + normalizedPosition[i] = + ( + ((baseMesh.points()[pointi]&userDirection) - startPosition) + / totalLength + ); + } - if (cosAngle < -0.9) - { - leftEdge[pointi] = edgei; - } - else if (cosAngle > 0.9) - { - rightEdge[pointi] = edgei; - } + // Sort the normalized position + labelList order; + sortedOrder(normalizedPosition, order); + + DynamicList seedPointDist; + + // Select points from finest refinement (one point-per plane) + scalar prevDist = GREAT; + forAll(order, i) + { + label pointi = order[i]; + scalar curDist = normalizedPosition[pointi]; + if (mag(curDist - prevDist) > meshRefiner_.mergeDistance()) + { + seedPointDist.append(curDist); + prevDist = curDist; } } + // Collect data from all processors + scalarList allSeedPointDist; + { + List gatheredDist(Pstream::nProcs()); + gatheredDist[Pstream::myProcNo()] = seedPointDist; + Pstream::gatherList(gatheredDist); - pointField newPoints(points); + // Combine processor lists into one big list. + allSeedPointDist = + ListListOps::combine + ( + gatheredDist, accessOp() + ); + } + + // Pre-set the points not to smooth (after expansion) + PackedBoolList isFrozenPoint(baseMesh.nPoints(), true); + + { + scalar minSeed = min(allSeedPointDist); + Pstream::scatter(minSeed); + scalar maxSeed = max(allSeedPointDist); + Pstream::scatter(maxSeed); + + forAll(normalizedPosition, posI) + { + const scalar pos = normalizedPosition[posI]; + if + ( + (mag(pos-minSeed) < meshRefiner_.mergeDistance()) + || (mag(pos-maxSeed) < meshRefiner_.mergeDistance()) + ) + { + // Boundary point: freeze + isFrozenPoint.set(pointLabels[posI]); + } + else + { + // Internal to moving region + isFrozenPoint.unset(pointLabels[posI]); + } + } + } Info<< "Smoothing " << geometry[surfi].name() << ':' << nl << " Direction : " << userDirection << nl << " Number of points : " << returnReduce(pointLabels.size(), sumOp