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; } } }