From 511431d6df33afcf3a53f396abc0a4445d726d76 Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 31 Jan 2022 16:43:00 +0000 Subject: [PATCH] BUG: snappyHexMesh: excessive memory blockLevel. Fixes #2345 Assumes that gap is formed when both surfaces agree i.e. it takes the minimum distance of the two. This means that any wave only needs to be propagated according to the originating surface. --- .../meshRefinement/meshRefinement.H | 27 -- .../meshRefinement/meshRefinementBaffles.C | 18 +- .../meshRefinement/meshRefinementBlock.C | 374 +++++++++--------- .../snappyHexMesh/meshRefinement/wallPoints.H | 15 +- .../meshRefinement/wallPointsI.H | 146 +++++-- 5 files changed, 322 insertions(+), 258 deletions(-) diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H index 0e753f206b..66fbf372fc 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H @@ -528,33 +528,6 @@ private: label& nRefine ) const; - void markMultiRegionCell - ( - const label celli, - //const label globalRegion, - //const label zonei, - const FixedList& surface, - - Map>& cellToRegions, - bitSet& isMultiRegion - ) const; - - void detectMultiRegionCells - ( - const labelListList& faceZones, - const labelList& testFaces, - - const labelList& surface1, - const List& hit1, - const labelList& region1, - - const labelList& surface2, - const List& hit2, - const labelList& region2, - - bitSet& isMultiRegion - ) const; - //- Mark cells for surface proximity based refinement. label markProximityRefinementWave ( diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C index 0f47b7e611..8bc18ac884 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBaffles.C @@ -6,7 +6,7 @@ \\/ M anipulation | ------------------------------------------------------------------------------- Copyright (C) 2011-2014 OpenFOAM Foundation - Copyright (C) 2015-2021 OpenCFD Ltd. + Copyright (C) 2015-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -57,6 +57,7 @@ License #include "FaceCellWave.H" #include "wallPoints.H" +#include "searchableSurfaces.H" // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // @@ -2409,8 +2410,21 @@ void Foam::meshRefinement::growCellZone List allFaceInfo(mesh_.nFaces()); + // No blocked faces, limitless gap size const bitSet isBlockedFace(mesh_.nFaces()); - wallPoints::trackData td(isBlockedFace); + List regionToBlockSize(surfaces_.surfaces().size()); + { + forAll(surfaces_.surfaces(), surfi) + { + const label geomi = surfaces_.surfaces()[surfi]; + const auto& s = surfaces_.geometry()[geomi]; + const label nRegions = s.regions().size(); + regionToBlockSize[surfi].setSize(nRegions, Foam::sqr(GREAT)); + } + } + + + wallPoints::trackData td(isBlockedFace, regionToBlockSize); FaceCellWave wallDistCalc ( mesh_, diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C index ff7e5ff556..7ddccb39da 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2018-2019 OpenCFD Ltd. + Copyright (C) 2018-2019,2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -294,151 +294,151 @@ void Foam::meshRefinement::growSet } -void Foam::meshRefinement::markMultiRegionCell -( - const label celli, - const FixedList& surface, - - Map>& cellToRegions, - bitSet& isMultiRegion -) const -{ - if (!isMultiRegion[celli]) - { - Map>::iterator fnd = cellToRegions.find(celli); - if (!fnd.found()) - { - cellToRegions.insert(celli, surface); - } - else if (fnd() != surface) - { - // Found multiple intersections on cell - isMultiRegion.set(celli); - } - } -} +//void Foam::meshRefinement::markMultiRegionCell +//( +// const label celli, +// const FixedList& surface, +// +// Map>& cellToRegions, +// bitSet& isMultiRegion +//) const +//{ +// if (!isMultiRegion[celli]) +// { +// Map>::iterator fnd = cellToRegions.find(celli); +// if (!fnd.found()) +// { +// cellToRegions.insert(celli, surface); +// } +// else if (fnd() != surface) +// { +// // Found multiple intersections on cell +// isMultiRegion.set(celli); +// } +// } +//} -void Foam::meshRefinement::detectMultiRegionCells -( - const labelListList& faceZones, - const labelList& testFaces, - - const labelList& surface1, - const List& hit1, - const labelList& region1, - - const labelList& surface2, - const List& hit2, - const labelList& region2, - - bitSet& isMultiRegion -) const -{ - isMultiRegion.clear(); - isMultiRegion.setSize(mesh_.nCells()); - - Map> cellToRegions(testFaces.size()); - - forAll(testFaces, i) - { - const pointIndexHit& info1 = hit1[i]; - if (info1.hit()) - { - const label facei = testFaces[i]; - const labelList& fz1 = faceZones[surface1[i]]; - const FixedList surfaceInfo1 - ({ - surface1[i], - region1[i], - (fz1.size() ? fz1[info1.index()] : region1[i]) - }); - - markMultiRegionCell - ( - mesh_.faceOwner()[facei], - surfaceInfo1, - cellToRegions, - isMultiRegion - ); - if (mesh_.isInternalFace(facei)) - { - markMultiRegionCell - ( - mesh_.faceNeighbour()[facei], - surfaceInfo1, - cellToRegions, - isMultiRegion - ); - } - - const pointIndexHit& info2 = hit2[i]; - - if (info2.hit() && info1 != info2) - { - const labelList& fz2 = faceZones[surface2[i]]; - const FixedList surfaceInfo2 - ({ - surface2[i], - region2[i], - (fz2.size() ? fz2[info2.index()] : region2[i]) - }); - - markMultiRegionCell - ( - mesh_.faceOwner()[facei], - surfaceInfo2, - cellToRegions, - isMultiRegion - ); - if (mesh_.isInternalFace(facei)) - { - markMultiRegionCell - ( - mesh_.faceNeighbour()[facei], - surfaceInfo2, - cellToRegions, - isMultiRegion - ); - } - } - } - } - - - if (debug&meshRefinement::MESH) - { - volScalarField multiCell - ( - IOobject - ( - "multiCell", - mesh_.time().timeName(), - mesh_, - IOobject::NO_READ, - IOobject::NO_WRITE, - false - ), - mesh_, - dimensionedScalar - ( - "zero", - dimensionSet(0, 1, 0, 0, 0), - 0.0 - ) - ); - forAll(isMultiRegion, celli) - { - if (isMultiRegion[celli]) - { - multiCell[celli] = 1.0; - } - } - - Info<< "Writing all multi cells to " << multiCell.name() << endl; - multiCell.write(); - } -} +//void Foam::meshRefinement::detectMultiRegionCells +//( +// const labelListList& faceZones, +// const labelList& testFaces, +// +// const labelList& surface1, +// const List& hit1, +// const labelList& region1, +// +// const labelList& surface2, +// const List& hit2, +// const labelList& region2, +// +// bitSet& isMultiRegion +//) const +//{ +// isMultiRegion.clear(); +// isMultiRegion.setSize(mesh_.nCells()); +// +// Map> cellToRegions(testFaces.size()); +// +// forAll(testFaces, i) +// { +// const pointIndexHit& info1 = hit1[i]; +// if (info1.hit()) +// { +// const label facei = testFaces[i]; +// const labelList& fz1 = faceZones[surface1[i]]; +// const FixedList surfaceInfo1 +// ({ +// surface1[i], +// region1[i], +// (fz1.size() ? fz1[info1.index()] : region1[i]) +// }); +// +// markMultiRegionCell +// ( +// mesh_.faceOwner()[facei], +// surfaceInfo1, +// cellToRegions, +// isMultiRegion +// ); +// if (mesh_.isInternalFace(facei)) +// { +// markMultiRegionCell +// ( +// mesh_.faceNeighbour()[facei], +// surfaceInfo1, +// cellToRegions, +// isMultiRegion +// ); +// } +// +// const pointIndexHit& info2 = hit2[i]; +// +// if (info2.hit() && info1 != info2) +// { +// const labelList& fz2 = faceZones[surface2[i]]; +// const FixedList surfaceInfo2 +// ({ +// surface2[i], +// region2[i], +// (fz2.size() ? fz2[info2.index()] : region2[i]) +// }); +// +// markMultiRegionCell +// ( +// mesh_.faceOwner()[facei], +// surfaceInfo2, +// cellToRegions, +// isMultiRegion +// ); +// if (mesh_.isInternalFace(facei)) +// { +// markMultiRegionCell +// ( +// mesh_.faceNeighbour()[facei], +// surfaceInfo2, +// cellToRegions, +// isMultiRegion +// ); +// } +// } +// } +// } +// +// +// if (debug&meshRefinement::MESH) +// { +// volScalarField multiCell +// ( +// IOobject +// ( +// "multiCell", +// mesh_.time().timeName(), +// mesh_, +// IOobject::NO_READ, +// IOobject::NO_WRITE, +// false +// ), +// mesh_, +// dimensionedScalar +// ( +// "zero", +// dimensionSet(0, 1, 0, 0, 0), +// 0.0 +// ) +// ); +// forAll(isMultiRegion, celli) +// { +// if (isMultiRegion[celli]) +// { +// multiCell[celli] = 1.0; +// } +// } +// +// Info<< "Writing all multi cells to " << multiCell.name() << endl; +// multiCell.write(); +// } +//} Foam::label Foam::meshRefinement::markProximityRefinementWave @@ -486,19 +486,18 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave // Re-work the blockLevel of the blockedSurfaces into a length scale // and a number of cells to cross - scalarField regionToBlockSize(surfaces_.blockLevel().size(), 0); - //label nIters = 2; - + List regionToBlockSize(surfaces_.surfaces().size()); for (const label surfi : blockedSurfaces) { const label geomi = surfaces_.surfaces()[surfi]; const searchableSurface& s = surfaces_.geometry()[geomi]; const label nRegions = s.regions().size(); + regionToBlockSize[surfi].setSize(nRegions); for (label regioni = 0; regioni < nRegions; regioni++) { const label globalRegioni = surfaces_.globalRegion(surfi, regioni); const label bLevel = surfaces_.blockLevel()[globalRegioni]; - regionToBlockSize[globalRegioni] = + regionToBlockSize[surfi][regioni] = meshCutter_.level0EdgeLength()/pow(2.0, bLevel); //const label mLevel = surfaces_.maxLevel()[globalRegioni]; @@ -513,6 +512,7 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave //nIters = max(nIters, (2<<(mLevel-bLevel))); } } + // Clever limiting of the number of iterations (= max cells in the channel) // since it has too many problematic issues, e.g. with volume refinement // and the real check uses the proper distance anyway just disable. @@ -540,7 +540,8 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave // TBD. correct nIters for actual cellLevel (since e.g. volume refinement // might add to cell level). Unfortunately we only have minLevel, // not maxLevel. Problem: what if volume refinement only in middle of - // channel? Workaround: have dummy surface with e.g. maxLevel 100 to + // channel? Say channel 1m wide with a 0.1m sphere of refinement + // Workaround: have dummy surface with e.g. maxLevel 100 to // force nIters to be high enough. @@ -577,22 +578,22 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave // Detect cells that are using multiple surface regions - bitSet isMultiRegion; - detectMultiRegionCells - ( - faceZones, - testFaces, - - surface1, - hit1, - region1, - - surface2, - hit2, - region2, - - isMultiRegion - ); + //bitSet isMultiRegion; + //detectMultiRegionCells + //( + // faceZones, + // testFaces, + // + // surface1, + // hit1, + // region1, + // + // surface2, + // hit2, + // region2, + // + // isMultiRegion + //); label n = 0; @@ -700,7 +701,7 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave // field. const bitSet isBlockedFace(intersectedFaces()); - wallPoints::trackData td(isBlockedFace); + wallPoints::trackData td(isBlockedFace, regionToBlockSize); FaceCellWave wallDistCalc ( mesh_, @@ -741,7 +742,7 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave { if (allCellInfo[celli].valid(wallDistCalc.data())) { - const point& cc = mesh_.C()[celli]; + const point& cc = mesh_.cellCentres()[celli]; // Nearest surface points const List& origin = allCellInfo[celli].origin(); @@ -792,7 +793,7 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave { if (allCellInfo[celli].valid(wallDistCalc.data())) { - const point& cc = mesh_.C()[celli]; + const point& cc = mesh_.cellCentres()[celli]; const List& origin = allCellInfo[celli].origin(); const List>& surface = @@ -803,15 +804,6 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave { for (label j = i + 1; j < origin.size(); j++) { - scalar maxDist - ( - max - ( - mag(cc-origin[i]), - mag(cc-origin[j]) - ) - ); - //if (isMultiRegion[celli]) //{ // // The intersection locations are too inaccurate @@ -830,21 +822,16 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave //else if (((cc-origin[i]) & (cc-origin[j])) < 0) { - const label globalRegioni = surfaces_.globalRegion - ( - surface[i][0], - surface[i][1] - ); - const label globalRegionj = surfaces_.globalRegion - ( - surface[j][0], - surface[j][1] - ); + const label surfi = surface[i][0]; + const label regioni = surface[i][1]; + + const label surfj = surface[j][0]; + const label regionj = surface[j][1]; const scalar maxSize = max ( - regionToBlockSize[globalRegioni], - regionToBlockSize[globalRegionj] + regionToBlockSize[surfi][regioni], + regionToBlockSize[surfj][regionj] ); if @@ -853,6 +840,15 @@ Foam::label Foam::meshRefinement::markProximityRefinementWave < Foam::sqr(2*maxSize) ) { + const scalar maxDist + ( + max + ( + mag(cc-origin[i]), + mag(cc-origin[j]) + ) + ); + smallGapDistance[celli] = max(smallGapDistance[celli], maxDist); nSmallGap++; diff --git a/src/mesh/snappyHexMesh/meshRefinement/wallPoints.H b/src/mesh/snappyHexMesh/meshRefinement/wallPoints.H index 015c5115de..aabb81a3c3 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/wallPoints.H +++ b/src/mesh/snappyHexMesh/meshRefinement/wallPoints.H @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2018-2020 OpenCFD Ltd. + Copyright (C) 2018-2020,2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -42,6 +42,7 @@ SourceFiles #include "tensor.H" #include "DynamicList.H" #include "labelList.H" +#include "scalarList.H" #include "bitSet.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -72,9 +73,17 @@ public: //- Per face whether the face should not be walked through const bitSet& isBlockedFace_; - trackData(const bitSet& isBlockedFace) + //- Per surface, per region the blockSize + const List& regionToBlockSize_; + + trackData + ( + const bitSet& isBlockedFace, + const List& regionToBlockSize + ) : - isBlockedFace_(isBlockedFace) + isBlockedFace_(isBlockedFace), + regionToBlockSize_(regionToBlockSize) {} }; diff --git a/src/mesh/snappyHexMesh/meshRefinement/wallPointsI.H b/src/mesh/snappyHexMesh/meshRefinement/wallPointsI.H index ddd9550367..4b2035131d 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/wallPointsI.H +++ b/src/mesh/snappyHexMesh/meshRefinement/wallPointsI.H @@ -5,7 +5,7 @@ \\ / A nd | www.openfoam.com \\/ M anipulation | ------------------------------------------------------------------------------- - Copyright (C) 2018-2020 OpenCFD Ltd. + Copyright (C) 2018-2022 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -202,26 +202,54 @@ inline bool Foam::wallPoints::updateCell const point& cc = mesh.cellCentres()[thisCelli]; bool hasChanged = false; + forAll(neighbourInfo.surface_, i) { const FixedList& nbrSurface = neighbourInfo.surface_[i]; - // Find in my surfaces - label index = surface_.find(nbrSurface); - if (index == -1) + // Check distance from nbr origin to cc against max walking distance + const scalar blockSize = + td.regionToBlockSize_[nbrSurface[0]][nbrSurface[1]]; + + const scalar d2 = magSqr(cc-neighbourInfo.origin_[i]); + + if (d2 < Foam::sqr(3*blockSize)) { - // Append - origin_.append(neighbourInfo.origin_[i]); - distSqr_.append(magSqr(cc-neighbourInfo.origin_[i])); - surface_.append(nbrSurface); - //normal_.append(neighbourInfo.normal_[i]); - hasChanged = true; + // Real distance less than max gap distance. Note that it should + // be at least 2 * blockSize (since gap is two cells across). + // Should be + // a little bit more to account for castellated path walk. Bit + // heuristic - try 3. Should be as low as possible to kill off + // unnecessary waves asap. + + // Find in my surfaces + label index = surface_.find(nbrSurface); + if (index == -1) + { + // Append + origin_.append(neighbourInfo.origin_[i]); + distSqr_.append(d2); + surface_.append(nbrSurface); + //normal_.append(neighbourInfo.normal_[i]); + hasChanged = true; + } + else + { + hasChanged = + update(cc, index, neighbourInfo, i, tol, td) + || hasChanged; + } } else { - hasChanged = - update(cc, index, neighbourInfo, i, tol, td) - || hasChanged; + // Real distance more than gap distance so ignore + //Pout<< "at cell:" << cc << " ignoring nbr info:" + // << neighbourInfo.origin_[i] + // << " from surface:" << nbrSurface[0] + // << " from region:" << nbrSurface[1] + // << " bloxkSize:" << blockSize + // << " distance:" << Foam::sqrt(d2) + // << endl; } } @@ -252,22 +280,44 @@ inline bool Foam::wallPoints::updateFace { const FixedList& nbrSurface = neighbourInfo.surface_[i]; - // Find in my surfaces - label index = surface_.find(nbrSurface); - if (index == -1) + // Check distance from nbr origin to cc against max walking distance + const scalar blockSize = + td.regionToBlockSize_[nbrSurface[0]][nbrSurface[1]]; + + const scalar d2 = magSqr(fc-neighbourInfo.origin_[i]); + + if (d2 < Foam::sqr(3*blockSize)) { - // Append - origin_.append(neighbourInfo.origin_[i]); - distSqr_.append(magSqr(fc-neighbourInfo.origin_[i])); - surface_.append(nbrSurface); - //normal_.append(neighbourInfo.normal_[i]); - hasChanged = true; + // Real distance less than max gap distance + + // Find in my surfaces + label index = surface_.find(nbrSurface); + if (index == -1) + { + // Append + origin_.append(neighbourInfo.origin_[i]); + distSqr_.append(d2); + surface_.append(nbrSurface); + //normal_.append(neighbourInfo.normal_[i]); + hasChanged = true; + } + else + { + hasChanged = + update(fc, index, neighbourInfo, i, tol, td) + || hasChanged; + } } else { - hasChanged = - update(fc, index, neighbourInfo, i, tol, td) - || hasChanged; + // Real distance more than gap distance so ignore + //Pout<< "at face:" << fc << " ignoring nbr info:" + // << neighbourInfo.origin_[i] + // << " from surface:" << nbrSurface[0] + // << " from region:" << nbrSurface[1] + // << " bloxkSize:" << blockSize + // << " distance:" << Foam::sqrt(d2) + // << endl; } } } @@ -298,22 +348,44 @@ inline bool Foam::wallPoints::updateFace { const FixedList& nbrSurface = neighbourInfo.surface_[i]; - // Find in my surfaces - const label index = surface_.find(nbrSurface); - if (index == -1) + // Check distance from nbr origin to cc against max walking distance + const scalar blockSize = + td.regionToBlockSize_[nbrSurface[0]][nbrSurface[1]]; + + const scalar d2 = magSqr(fc-neighbourInfo.origin_[i]); + + if (d2 < Foam::sqr(3*blockSize)) { - // Append - origin_.append(neighbourInfo.origin_[i]); - distSqr_.append(magSqr(fc-neighbourInfo.origin_[i])); - surface_.append(nbrSurface); - //normal_.append(neighbourInfo.normal_[i]); - hasChanged = true; + // Real distance less than max gap distance + + // Find in my surfaces + const label index = surface_.find(nbrSurface); + if (index == -1) + { + // Append + origin_.append(neighbourInfo.origin_[i]); + distSqr_.append(d2); + surface_.append(nbrSurface); + //normal_.append(neighbourInfo.normal_[i]); + hasChanged = true; + } + else + { + hasChanged = + update(fc, index, neighbourInfo, i, tol, td) + || hasChanged; + } } else { - hasChanged = - update(fc, index, neighbourInfo, i, tol, td) - || hasChanged; + // Real distance more than gap distance so ignore + //Pout<< "at face:" << fc << " ignoring nbr info:" + // << neighbourInfo.origin_[i] + // << " from surface:" << nbrSurface[0] + // << " from region:" << nbrSurface[1] + // << " bloxkSize:" << blockSize + // << " distance:" << Foam::sqrt(d2) + // << endl; } } }