From b7c54bc00c82c5a136599fb94c027c972a6dd57c Mon Sep 17 00:00:00 2001 From: mattijs Date: Mon, 9 Dec 2019 14:29:21 +0000 Subject: [PATCH] ENH: snappyHexMesh: proximity check This adds automatic deletion of cells inside small gaps. This is generally used to avoid having excessive numbers of cells in irrelevant areas of a geometry. It is nearly the opposite of automatic gap refinement - that refines cells to resolve the gap; this functionality removes cells to not mesh the gap. The proximity handling will remove those cells which are inside 'thin' gaps where 'thin' is defined as a distance of 2*'blockLevel' It will - detect surfaces which have the new 'blockLevel' specification - convert this to a minimum gap distance - detect cells which are inside this gap - remove these cells and add exposed faces to the nearest 'real' patch --- src/mesh/snappyHexMesh/Make/files | 1 + src/mesh/snappyHexMesh/Make/options | 3 +- .../meshRefinement/meshRefinement.H | 41 ++ .../meshRefinement/meshRefinementBlock.C | 685 +++++++++++++++++- .../snappyHexMesh/meshRefinement/wallPoints.C | 56 ++ .../snappyHexMesh/meshRefinement/wallPoints.H | 260 +++++++ .../meshRefinement/wallPointsI.H | 352 +++++++++ .../snappyHexMeshDriver/snappyRefineDriver.C | 25 + 8 files changed, 1418 insertions(+), 5 deletions(-) create mode 100644 src/mesh/snappyHexMesh/meshRefinement/wallPoints.C create mode 100644 src/mesh/snappyHexMesh/meshRefinement/wallPoints.H create mode 100644 src/mesh/snappyHexMesh/meshRefinement/wallPointsI.H diff --git a/src/mesh/snappyHexMesh/Make/files b/src/mesh/snappyHexMesh/Make/files index 0bdadb9848..2b5f896ec4 100644 --- a/src/mesh/snappyHexMesh/Make/files +++ b/src/mesh/snappyHexMesh/Make/files @@ -15,6 +15,7 @@ meshRefinement/meshRefinementProblemCells.C meshRefinement/meshRefinementRefine.C meshRefinement/meshRefinementGapRefine.C meshRefinement/meshRefinementBlock.C +meshRefinement/wallPoints.C meshRefinement/patchFaceOrientation.C refinementFeatures/refinementFeatures.C diff --git a/src/mesh/snappyHexMesh/Make/options b/src/mesh/snappyHexMesh/Make/options index 519e861b82..a5fb38a59e 100644 --- a/src/mesh/snappyHexMesh/Make/options +++ b/src/mesh/snappyHexMesh/Make/options @@ -7,7 +7,8 @@ EXE_INC = \ -I$(LIB_SRC)/sampling/lnInclude \ -I$(LIB_SRC)/lagrangian/basic/lnInclude \ -I$(LIB_SRC)/overset/lnInclude \ - -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude + -I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \ + -I$(LIB_SRC)/parallel/distributed/lnInclude LIB_LIBS = \ -lfiniteVolume \ diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H index 41766f7f20..3a43723a5e 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H @@ -533,6 +533,47 @@ 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 + ( + const scalar planarCos, + const labelList& blockedSurfaces, + const label nAllowRefine, + const labelList& neiLevel, + const pointField& neiCc, + + labelList& refineCell, + label& nRefine + ) const; + + // Baffle handling //- Get faces to repatch. Returns map from face to patch. diff --git a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C index 8ae862a1aa..2eb37d76f8 100644 --- a/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C +++ b/src/mesh/snappyHexMesh/meshRefinement/meshRefinementBlock.C @@ -33,6 +33,12 @@ License #include "unitConversion.H" #include "bitSet.H" +#include "FaceCellWave.H" +#include "volFields.H" +#include "wallPoints.H" +#include "searchableSurfaces.H" +#include "distributedTriSurfaceMesh.H" + // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // //Foam::label Foam::meshRefinement::markFakeGapRefinement @@ -288,6 +294,636 @@ 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::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 +( + const scalar planarCos, + const labelList& blockedSurfaces, + const label nAllowRefine, + const labelList& neiLevel, + const pointField& neiCc, + + labelList& refineCell, + label& nRefine +) const +{ + labelListList faceZones(surfaces_.surfaces().size()); + { + // If triSurface do additional zoning based on connectivity + for (const label surfi : blockedSurfaces) + { + const label geomi = surfaces_.surfaces()[surfi]; + const searchableSurface& s = surfaces_.geometry()[geomi]; + if (isA(s) && !isA(s)) + { + const triSurfaceMesh& surf = refCast(s); + const labelListList& edFaces = surf.edgeFaces(); + boolList isOpenEdge(edFaces.size(), false); + forAll(edFaces, i) + { + if (edFaces[i].size() == 1) + { + isOpenEdge[i] = true; + } + } + + labelList faceZone; + const label nZones = surf.markZones(isOpenEdge, faceZone); + if (nZones > 1) + { + faceZones[surfi].transfer(faceZone); + } + } + } + } + + + // 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; + + for (const label surfi : blockedSurfaces) + { + const label geomi = surfaces_.surfaces()[surfi]; + const searchableSurface& s = surfaces_.geometry()[geomi]; + const label nRegions = s.regions().size(); + for (label regioni = 0; regioni < nRegions; regioni++) + { + const label globalRegioni = surfaces_.globalRegion(surfi, regioni); + const label bLevel = surfaces_.blockLevel()[globalRegioni]; + regionToBlockSize[globalRegioni] = + meshCutter_.level0EdgeLength()/pow(2.0, bLevel); + + //const label mLevel = surfaces_.maxLevel()[globalRegioni]; + //// TBD: check for higher cached level of surface due to vol + //// refinement. Problem: might still miss refinement bubble + //// fully inside thin channel + //if (isA(s) && !isA(s)) + //{ + // const triSurfaceMesh& surf = refCast(s); + //} + + //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. + // Since the real check uses the proper distance anyway just disable. + const label nIters = mesh_.globalData().nTotalCells(); + + + // Collect candidate faces (i.e. intersecting any surface and + // owner/neighbour not yet refined) + const labelList testFaces(getRefineCandidateFaces(refineCell)); + + // Collect segments + pointField start(testFaces.size()); + pointField end(testFaces.size()); + labelList minLevel(testFaces.size()); + + calcCellCellRays + ( + neiCc, + neiLevel, + testFaces, + start, + end, + minLevel + ); + // 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 + // force nIters to be high enough. + + + // Test for all intersections (with surfaces of higher gap level than + // minLevel) and cache per cell the max surface level and the local normal + // on that surface. + + labelList surface1; + List hit1; + labelList region1; + vectorField normal1; + + labelList surface2; + List hit2; + labelList region2; + vectorField normal2; + + surfaces_.findNearestIntersection + ( + blockedSurfaces, + start, + end, + + surface1, + hit1, + region1, // local region + normal1, + + surface2, + hit2, + region2, // local region + normal2 + ); + + + // Detect cells that are using multiple surface regions + bitSet isMultiRegion; + detectMultiRegionCells + ( + faceZones, + testFaces, + + surface1, + hit1, + region1, + + surface2, + hit2, + region2, + + isMultiRegion + ); + + + label n = 0; + forAll(testFaces, i) + { + if (hit1[i].hit()) + { + n++; + } + } + + List faceDist(n); + labelList changedFaces(n); + n = 0; + + DynamicList originLocation(2); + DynamicList originDistSqr(2); + DynamicList> originSurface(2); + //DynamicList originNormal(2); + + forAll(testFaces, i) + { + if (hit1[i].hit()) + { + const label facei = testFaces[i]; + const point& fc = mesh_.faceCentres()[facei]; + const labelList& fz1 = faceZones[surface1[i]]; + + originLocation.clear(); + originDistSqr.clear(); + //originNormal.clear(); + originSurface.clear(); + + originLocation.append(hit1[i].hitPoint()); + originDistSqr.append(magSqr(fc-hit1[i].hitPoint())); + //originNormal.append(normal1[i]); + originSurface.append + ( + FixedList + ({ + surface1[i], + region1[i], + (fz1.size() ? fz1[hit1[i].index()] : region1[i]) + }) + ); + + if (hit2[i].hit() && hit1[i] != hit2[i]) + { + const labelList& fz2 = faceZones[surface2[i]]; + originLocation.append(hit2[i].hitPoint()); + originDistSqr.append(magSqr(fc-hit2[i].hitPoint())); + //originNormal.append(normal2[i]); + originSurface.append + ( + FixedList + ({ + surface2[i], + region2[i], + (fz2.size() ? fz2[hit2[i].index()] : region2[i]) + }) + ); + } + + faceDist[n] = wallPoints + ( + originLocation, // origin + originDistSqr, // distance to origin + originSurface // surface+region+zone + //originNormal // normal at origin + ); + changedFaces[n] = facei; + n++; + } + } + + + // Clear intersection info + surface1.clear(); + hit1.clear(); + region1.clear(); + normal1.clear(); + surface2.clear(); + hit2.clear(); + region2.clear(); + normal2.clear(); + + + List allFaceInfo(mesh_.nFaces()); + List allCellInfo(mesh_.nCells()); + + FaceCellWave wallDistCalc + ( + mesh_, + changedFaces, + faceDist, + allFaceInfo, + allCellInfo, + 0 // max iterations + ); + wallDistCalc.iterate(nIters); + + + if (debug&meshRefinement::MESH) + { + // Dump current nearest opposite surfaces + volScalarField distance + ( + IOobject + ( + "gapSize", + mesh_.time().timeName(), + mesh_, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh_, + dimensionedScalar + ( + "zero", + dimLength, //dimensionSet(0, 1, 0, 0, 0), + -1 + ) + ); + + forAll(allCellInfo, celli) + { + if (allCellInfo[celli].valid(wallDistCalc.data())) + { + const point& cc = mesh_.C()[celli]; + // Nearest surface points + const List& origin = allCellInfo[celli].origin(); + + // Find 'opposite' pair with minimum distance + for (label i = 0; i < origin.size(); i++) + { + for (label j = i + 1; j < origin.size(); j++) + { + if (((cc-origin[i]) & (cc-origin[j])) < 0) + { + const scalar d(mag(origin[i]-origin[j])); + if (distance[celli] < 0) + { + distance[celli] = d; + } + else + { + distance[celli] = min(distance[celli], d); + } + } + } + } + } + } + distance.correctBoundaryConditions(); + + Info<< "Writing measured gap distance to " + << distance.name() << endl; + distance.write(); + } + + + + // Detect tight gaps: + // - cell is inbetween the two surfaces + // - two surfaces are planarish + // - two surfaces are not too far apart + // (number of walking iterations is a too-coarse measure) + + scalarField smallGapDistance(mesh_.nCells(), 0.0); + label nMulti = 0; + label nSmallGap = 0; + + //OBJstream str(mesh_.time().timePath()/"multiRegion.obj"); + + + forAll(allCellInfo, celli) + { + if (allCellInfo[celli].valid(wallDistCalc.data())) + { + const point& cc = mesh_.C()[celli]; + + const List& origin = allCellInfo[celli].origin(); + const List>& surface = + allCellInfo[celli].surface(); + + // Find pair with minimum distance + for (label i = 0; i < origin.size(); i++) + { + 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 + // // (since not proper nearest, just a cell-cell ray + // // intersection) so include always + // + // smallGapDistance[celli] = + // max(smallGapDistance[celli], maxDist); + // + // + // str.write(linePointRef(cc, origin[i])); + // str.write(linePointRef(cc, origin[j])); + // + // nMulti++; + //} + //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 scalar maxSize = max + ( + regionToBlockSize[globalRegioni], + regionToBlockSize[globalRegionj] + ); + + if + ( + magSqr(origin[i]-origin[j]) + < Foam::sqr(2*maxSize) + ) + { + smallGapDistance[celli] = + max(smallGapDistance[celli], maxDist); + nSmallGap++; + } + } + } + } + } + } + + + if (debug) + { + Info<< "Marked for blocking due to intersecting multiple surfaces : " + << returnReduce(nMulti, sumOp