zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

snappyHexMesh : refine based on curvature

Browse Source
This commit is contained in:
Mattijs Janssens
2022-08-04 17:09:38 +00:00
parent 227727d413
commit 27c3d0c23b
20 changed files with 2326 additions and 35 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
src/mesh/snappyHexMesh/meshRefinement/meshRefinement.H
View File

@ -471,6 +471,19 @@ private:
label& nRefine
) const;
//- Refine cells containing triangles with high refinement level
// (currently due to high curvature or being inside shell with
// high level)
label markSurfaceFieldRefinement
(
const label nAllowRefine,
const labelList& neiLevel,
const pointField& neiCc,
labelList& refineCell,
label& nRefine
) const;
//- Helper: count number of normals1 that are in normals2
label countMatches
(
@ -479,6 +492,18 @@ private:
const scalar tol = 1e-6
) const;
//- Detect if two intersection points are high curvature (w.r.t.
// lengthscale
bool highCurvature
(
const scalar minCosAngle,
const scalar lengthScale,
const point& p0,
const vector& n0,
const point& p1,
const vector& n1
) const;
//- Mark cells for surface curvature based refinement. Marks if
// local curvature > curvature or if on different regions
// (markDifferingRegions)

116
src/mesh/snappyHexMesh/meshRefinement/meshRefinementGapRefine.C
View File

@ -1833,4 +1833,120 @@ Foam::label Foam::meshRefinement::markSmallFeatureRefinement
}
////XXXXXXXX
Foam::label Foam::meshRefinement::markSurfaceFieldRefinement
(
const label nAllowRefine,
const labelList& neiLevel,
const pointField& neiCc,
labelList& refineCell,
label& nRefine
) const
{
const labelList& cellLevel = meshCutter_.cellLevel();
const labelList& surfaceIndices = surfaces_.surfaces();
label oldNRefine = nRefine;
//- Force calculation of tetBasePt
(void)mesh_.tetBasePtIs();
(void)mesh_.cellTree();
const indexedOctree<treeDataCell>& tree = mesh_.cellTree();
forAll(surfaceIndices, surfI)
{
label geomI = surfaceIndices[surfI];
const searchableSurface& geom = surfaces_.geometry()[geomI];
// Get the element index in a roundabout way. Problem is e.g.
// distributed surface where local indices differ from global
// ones (needed for getRegion call)
pointField ctrs;
labelList region;
labelList minLevelField;
{
// Representative local coordinates and bounding sphere
scalarField radiusSqr;
geom.boundingSpheres(ctrs, radiusSqr);
List<pointIndexHit> info;
geom.findNearest(ctrs, radiusSqr, info);
forAll(info, i)
{
if (!info[i].hit())
{
FatalErrorInFunction
<< "fc:" << ctrs[i]
<< " radius:" << radiusSqr[i]
<< exit(FatalError);
}
}
geom.getRegion(info, region);
geom.getField(info, minLevelField);
}
if (minLevelField.size() != geom.size())
{
Pout<< "** no minLevelField" << endl;
continue;
}
label nOldRefine = 0;
forAll(ctrs, i)
{
label cellI = -1;
if (tree.nodes().size() && tree.bb().contains(ctrs[i]))
{
cellI = tree.findInside(ctrs[i]);
}
if
(
cellI != -1
&& refineCell[cellI] == -1
&& minLevelField[i] > cellLevel[cellI]
)
{
if
(
!markForRefine
(
surfI,
nAllowRefine,
refineCell[cellI],
nRefine
)
)
{
break;
}
}
}
Info<< "For surface " << geom.name() << " found "
<< returnReduce(nRefine-nOldRefine, sumOp<label>())
<< " cells containing cached refinement field" << endl;
if
(
returnReduce(nRefine, sumOp<label>())
> returnReduce(nAllowRefine, sumOp<label>())
)
{
Info<< "Reached refinement limit." << endl;
}
}
return returnReduce(nRefine-oldNRefine, sumOp<label>());
}
////XXXXXXXXX
// ************************************************************************* //

202
src/mesh/snappyHexMesh/meshRefinement/meshRefinementRefine.C
View File

@ -1109,6 +1109,71 @@ Foam::label Foam::meshRefinement::countMatches
}
//XXXXXX
//bool Foam::meshRefinement::highCurvature
//(
// const scalar minCosAngle,
// const point& p0,
// const vector& n0,
// const point& p1,
// const vector& n1
//) const
//{
// return ((n0&n1) < minCosAngle);
//}
bool Foam::meshRefinement::highCurvature
(
const scalar minCosAngle,
const scalar lengthScale,
const point& p0,
const vector& n0,
const point& p1,
const vector& n1
) const
{
const scalar cosAngle = (n0&n1);
if (cosAngle < minCosAngle)
{
// Sharp feature, independent of intersection points
return true;
}
else if (cosAngle > 1-1e-6)
{
// Co-planar
return false;
}
else
{
// Calculate radius of curvature
vector axis(n0 ^ n1);
const plane pl0(p0, (n0 ^ axis));
const scalar r1 = pl0.normalIntersect(p1, n1);
//const point radiusPoint(p1+r1*n1);
//DebugVar(radiusPoint);
const plane pl1(p1, (n1 ^ axis));
const scalar r0 = pl1.normalIntersect(p0, n0);
//const point radiusPoint(p0+r0*n0);
//DebugVar(radiusPoint);
//- Take average radius. Bit ad hoc
const scalar radius = 0.5*(mag(r1)+mag(r0));
if (radius < lengthScale)
{
return true;
}
else
{
return false;
}
}
}
//XXXXX
// Mark cells for surface curvature based refinement.
Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
(
@ -1187,9 +1252,12 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
// Test for all intersections (with surfaces of higher max level than
// minLevel) and cache per cell the interesting inter
labelListList cellSurfLevels(mesh_.nCells());
List<pointList> cellSurfLocations(mesh_.nCells());
List<vectorList> cellSurfNormals(mesh_.nCells());
{
// Per segment the intersection point of the surfaces
List<pointList> surfaceLocation;
// Per segment the normals of the surfaces
List<vectorList> surfaceNormal;
// Per segment the list of levels of the surfaces
@ -1205,6 +1273,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
labelList(surfaces_.maxLevel().size(), 0), // min level
surfaces_.maxLevel(),
surfaceLocation,
surfaceNormal,
surfaceLevel
);
@ -1216,12 +1285,14 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
labelList visitOrder;
forAll(surfaceNormal, pointI)
{
pointList& pLocations = surfaceLocation[pointI];
vectorList& pNormals = surfaceNormal[pointI];
labelList& pLevel = surfaceLevel[pointI];
sortedOrder(pNormals, visitOrder, normalLess(pNormals));
sortedOrder(pNormals, visitOrder, normalLess(pLocations));
pNormals = List<point>(pNormals, visitOrder);
pLocations = List<point>(pLocations, visitOrder);
pNormals = List<vector>(pNormals, visitOrder);
pLevel = labelUIndList(pLevel, visitOrder);
}
@ -1236,6 +1307,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
label faceI = testFaces[i];
label own = mesh_.faceOwner()[faceI];
const pointList& fPoints = surfaceLocation[i];
const vectorList& fNormals = surfaceNormal[i];
const labelList& fLevels = surfaceLevel[i];
@ -1244,6 +1316,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
if (fLevels[hitI] > cellLevel[own])
{
cellSurfLevels[own].append(fLevels[hitI]);
cellSurfLocations[own].append(fPoints[hitI]);
cellSurfNormals[own].append(fNormals[hitI]);
}
@ -1253,6 +1326,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
if (fLevels[hitI] > cellLevel[nei])
{
cellSurfLevels[nei].append(fLevels[hitI]);
cellSurfLocations[nei].append(fPoints[hitI]);
cellSurfNormals[nei].append(fNormals[hitI]);
}
}
@ -1266,7 +1340,7 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
if (debug)
{
label nSet = 0;
label nNormals = 9;
label nNormals = 0;
forAll(cellSurfNormals, cellI)
{
const vectorList& normals = cellSurfNormals[cellI];
@ -1296,21 +1370,38 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
for
(
label cellI = 0;
!reachedLimit && cellI < cellSurfNormals.size();
!reachedLimit && cellI < cellSurfLocations.size();
cellI++
)
{
const pointList& points = cellSurfLocations[cellI];
const vectorList& normals = cellSurfNormals[cellI];
const labelList& levels = cellSurfLevels[cellI];
// Current cell size
const scalar cellSize =
meshCutter_.level0EdgeLength()/pow(2.0, cellLevel[cellI]);
// n^2 comparison of all normals in a cell
for (label i = 0; !reachedLimit && i < normals.size(); i++)
{
for (label j = i+1; !reachedLimit && j < normals.size(); j++)
{
if ((normals[i] & normals[j]) < curvature)
//if ((normals[i] & normals[j]) < curvature)
if
(
highCurvature
(
curvature,
cellSize,
points[i],
normals[i],
points[j],
normals[j]
)
)
{
label maxLevel = max(levels[i], levels[j]);
const label maxLevel = max(levels[i], levels[j]);
if (cellLevel[cellI] < maxLevel)
{
@ -1358,8 +1449,10 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
label own = mesh_.faceOwner()[faceI];
label nei = mesh_.faceNeighbour()[faceI];
const pointList& ownPoints = cellSurfLocations[own];
const vectorList& ownNormals = cellSurfNormals[own];
const labelList& ownLevels = cellSurfLevels[own];
const pointList& neiPoints = cellSurfLocations[nei];
const vectorList& neiNormals = cellSurfNormals[nei];
const labelList& neiLevels = cellSurfLevels[nei];
@ -1379,13 +1472,30 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
if (!ownIsSubset && !neiIsSubset)
{
// Current average cell size. Is min? max? average?
const scalar cellSize =
meshCutter_.level0EdgeLength()
/ pow(2.0, min(cellLevel[own], cellLevel[nei]));
// n^2 comparison of between ownNormals and neiNormals
for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
{
for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
{
// Have valid data on both sides. Check curvature.
if ((ownNormals[i] & neiNormals[j]) < curvature)
//if ((ownNormals[i] & neiNormals[j]) < curvature)
if
(
highCurvature
(
curvature,
cellSize,
ownPoints[i],
ownNormals[i],
neiPoints[j],
neiNormals[j]
)
)
{
// See which side to refine.
if (cellLevel[own] < ownLevels[i])
@ -1441,7 +1551,11 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
}
// Send over surface normal to neighbour cell.
// Send over surface point/normal to neighbour cell.
// labelListList neiSurfaceLevel;
// syncTools::swapBoundaryCellList(mesh_, cellSurfLevels, neiSurfaceLevel);
List<vectorList> neiSurfacePoints;
syncTools::swapBoundaryCellList(mesh_, cellSurfLocations, neiSurfacePoints);
List<vectorList> neiSurfaceNormals;
syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
@ -1456,9 +1570,13 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
label own = mesh_.faceOwner()[faceI];
label bFaceI = faceI - mesh_.nInternalFaces();
const pointList& ownPoints = cellSurfLocations[own];
const vectorList& ownNormals = cellSurfNormals[own];
const labelList& ownLevels = cellSurfLevels[own];
const pointList& neiPoints = neiSurfacePoints[bFaceI];
const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
//const labelList& neiLevels = neiSurfaceLevel[bFacei];
// Special case: owner normals set is a subset of the neighbour
// normals. Do not do curvature refinement since other cell's normals
@ -1475,13 +1593,30 @@ Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
if (!ownIsSubset && !neiIsSubset)
{
// Current average cell size. Is min? max? average?
const scalar cellSize =
meshCutter_.level0EdgeLength()
/ pow(2.0, min(cellLevel[own], neiLevel[bFaceI]));
// n^2 comparison of between ownNormals and neiNormals
for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
{
for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
{
// Have valid data on both sides. Check curvature.
if ((ownNormals[i] & neiNormals[j]) < curvature)
//if ((ownNormals[i] & neiNormals[j]) < curvature)
if
(
highCurvature
(
curvature,
cellSize,
ownPoints[i],
ownNormals[i],
neiPoints[j],
neiNormals[j]
)
)
{
if (cellLevel[own] < ownLevels[i])
{
@ -2205,25 +2340,46 @@ Foam::labelList Foam::meshRefinement::refineCandidates
// Refinement based on features smaller than cell size
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if
(
smallFeatureRefinement
&& (planarCos >= -1 && planarCos <= 1)
&& max(shells_.maxGapLevel()) > 0
)
if (smallFeatureRefinement)
{
label nGap = markSmallFeatureRefinement
label nGap = 0;
if
(
planarCos,
nAllowRefine,
neiLevel,
neiCc,
planarCos >= -1
&& planarCos <= 1
&& max(shells_.maxGapLevel()) > 0
)
{
nGap = markSmallFeatureRefinement
(
planarCos,
nAllowRefine,
neiLevel,
neiCc,
refineCell,
nRefine
);
refineCell,
nRefine
);
}
Info<< "Marked for refinement due to close opposite surfaces "
<< ": " << nGap << " cells." << endl;
label nCurv = 0;
if (max(surfaces_.maxCurvatureLevel()) > 0)
{
nCurv = markSurfaceFieldRefinement
(
nAllowRefine,
neiLevel,
neiCc,
refineCell,
nRefine
);
Info<< "Marked for refinement"
<< " due to curvature "
<< ": " << nCurv << " cells." << endl;
}
}

379
src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
Copyright (C) 2015-2020 OpenCFD Ltd.
Copyright (C) 2015-2022 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -38,6 +38,8 @@ License
// For dictionary::get wrapper
#include "meshRefinement.H"
#include "OBJstream.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::labelList Foam::refinementSurfaces::findHigherLevel
@ -194,15 +196,14 @@ Foam::refinementSurfaces::refinementSurfaces
labelList globalMaxLevel(surfI, Zero);
labelList globalLevelIncr(surfI, Zero);
FixedList<label, 3> nullGapLevel;
nullGapLevel[0] = 0;
nullGapLevel[1] = 0;
nullGapLevel[2] = 0;
const FixedList<label, 3> nullGapLevel({0, 0, 0});
List<FixedList<label, 3>> globalGapLevel(surfI);
List<volumeType> globalGapMode(surfI);
boolList globalGapSelf(surfI);
const FixedList<label, 4> nullCurvLevel({0, 0, 0, -1});
List<FixedList<label, 4>> globalCurvLevel(surfI);
scalarField globalAngle(surfI, -GREAT);
PtrList<dictionary> globalPatchInfo(surfI);
@ -216,6 +217,7 @@ Foam::refinementSurfaces::refinementSurfaces
List<Map<FixedList<label, 3>>> regionGapLevel(surfI);
List<Map<volumeType>> regionGapMode(surfI);
List<Map<bool>> regionGapSelf(surfI);
List<Map<FixedList<label, 4>>> regionCurvLevel(surfI);
List<Map<scalar>> regionAngle(surfI);
List<Map<autoPtr<dictionary>>> regionPatchInfo(surfI);
List<Map<label>> regionBlockLevel(surfI);
@ -307,6 +309,19 @@ Foam::refinementSurfaces::refinementSurfaces
globalGapSelf[surfI] =
dict.getOrDefault<bool>("gapSelf", true);
globalCurvLevel[surfI] = nullCurvLevel;
if (dict.readIfPresent("curvatureLevel", globalCurvLevel[surfI]))
{
if (globalCurvLevel[surfI][0] <= 0)
{
FatalIOErrorInFunction(dict)
<< "Illegal curvatureLevel specification for surface "
<< names_[surfI]
<< " : curvatureLevel:" << globalCurvLevel[surfI]
<< exit(FatalIOError);
}
}
const searchableSurface& surface = allGeometry_[surfaces_[surfI]];
// Surface zones
@ -435,6 +450,20 @@ Foam::refinementSurfaces::refinementSurfaces
)
);
FixedList<label, 4> curvSpec(nullCurvLevel);
if (regionDict.readIfPresent("curvatureLevel", curvSpec))
{
if (curvSpec[0] <= 0)
{
FatalIOErrorInFunction(dict)
<< "Illegal curvatureLevel specification for surface "
<< names_[surfI]
<< " : curvatureLevel:" << curvSpec
<< exit(FatalIOError);
}
}
regionCurvLevel[surfI].insert(regionI, curvSpec);
if (regionDict.found("perpendicularAngle"))
{
regionAngle[surfI].insert
@ -505,6 +534,8 @@ Foam::refinementSurfaces::refinementSurfaces
extendedGapMode_ = volumeType::UNKNOWN;
selfProximity_.setSize(nRegions);
selfProximity_ = true;
extendedCurvatureLevel_.setSize(nRegions);
extendedCurvatureLevel_ = nullCurvLevel;
perpendicularAngle_.setSize(nRegions);
perpendicularAngle_ = -GREAT;
patchInfo_.setSize(nRegions);
@ -531,6 +562,8 @@ Foam::refinementSurfaces::refinementSurfaces
extendedGapLevel_[globalRegionI] = globalGapLevel[surfI];
extendedGapMode_[globalRegionI] = globalGapMode[surfI];
selfProximity_[globalRegionI] = globalGapSelf[surfI];
extendedCurvatureLevel_[globalRegionI] =
globalCurvLevel[surfI];
perpendicularAngle_[globalRegionI] = globalAngle[surfI];
if (globalPatchInfo.set(surfI))
{
@ -560,6 +593,8 @@ Foam::refinementSurfaces::refinementSurfaces
regionGapMode[surfI][iter.key()];
selfProximity_[globalRegionI] =
regionGapSelf[surfI][iter.key()];
extendedCurvatureLevel_[globalRegionI] =
regionCurvLevel[surfI][iter.key()];
}
forAllConstIters(regionAngle[surfI], iter)
{
@ -713,6 +748,26 @@ Foam::labelList Foam::refinementSurfaces::maxGapLevel() const
}
Foam::labelList Foam::refinementSurfaces::maxCurvatureLevel() const
{
labelList surfaceMax(surfaces_.size(), Zero);
forAll(surfaces_, surfI)
{
const wordList& regionNames = allGeometry_[surfaces_[surfI]].regions();
forAll(regionNames, regionI)
{
label globalI = globalRegion(surfI, regionI);
const FixedList<label, 4>& gapInfo =
extendedCurvatureLevel_[globalI];
surfaceMax[surfI] = max(surfaceMax[surfI], gapInfo[2]);
}
}
return surfaceMax;
}
// Precalculate the refinement level for every element of the searchable
// surface.
void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
@ -763,13 +818,15 @@ void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
// In case of triangulated surfaces only cache value if triangle
// centre and vertices are in same shell
if (isA<triSurface>(geom))
const auto* tsPtr = isA<const triSurface>(geom);
if (tsPtr)
{
label nUncached = 0;
// Check if points differing from ctr level
const triSurface& ts = refCast<const triSurface>(geom);
const triSurface& ts = *tsPtr;
const pointField& points = ts.points();
// Determine minimum expected level to avoid having to
@ -845,6 +902,312 @@ void Foam::refinementSurfaces::setMinLevelFields(const shellSurfaces& shells)
}
// Precalculate the refinement level for every element of the searchable
// surface.
void Foam::refinementSurfaces::setCurvatureMinLevelFields
(
const scalar cosAngle,
const scalar level0EdgeLength
)
{
const labelList maxCurvLevel(maxCurvatureLevel());
forAll(surfaces_, surfI)
{
// Check if there is a specification of the extended curvature
if (maxCurvLevel[surfI] <= 0)
{
continue;
}
const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
const auto* tsPtr = isA<const triSurface>(geom);
// Cache the refinement level (max of surface level and shell level)
// on a per-element basis. Only makes sense if there are lots of
// elements. Possibly should have 'enough' elements to have fine
// enough resolution but for now just make sure we don't catch e.g.
// searchableBox (size=6)
if (tsPtr)
{
// Representative local coordinates and bounding sphere
pointField ctrs;
scalarField radiusSqr;
geom.boundingSpheres(ctrs, radiusSqr);
labelList minLevelField(ctrs.size(), Zero);
//labelList surfMin(ctrs.size(), Zero);
labelList surfMax(ctrs.size(), Zero);
labelList nCurvCells(ctrs.size(), Zero);
labelList curvIgnore(ctrs.size(), -1);
{
// Get the element index in a roundabout way. Problem is e.g.
// distributed surface where local indices differ from global
// ones (needed for getRegion call)
List<pointIndexHit> info;
geom.findNearest(ctrs, radiusSqr, info);
// Get per element the region
labelList region;
geom.getRegion(info, region);
// See if a cached level field available (from e.g. shells)
labelList cachedField;
geom.getField(info, cachedField);
// From the region get the surface-wise refinement level
forAll(minLevelField, i)
{
if (info[i].hit()) //Note: should not be necessary
{
const label globali = globalRegion(surfI, region[i]);
curvIgnore[i] = extendedCurvatureLevel_[globali][3];
nCurvCells[i] = extendedCurvatureLevel_[globali][0];
//surfMin[i] = extendedCurvatureLevel_[globali][1];
surfMax[i] = extendedCurvatureLevel_[globali][2];
minLevelField[i] = minLevel(surfI, region[i]);
if (cachedField.size() > i)
{
minLevelField[i] =
max(minLevelField[i], cachedField[i]);
}
}
}
}
// Calculate per-triangle curvature. This is the max of the
// measured point-based curvature + some constraints.
scalarField cellCurv(ctrs.size(), Zero);
{
// Walk surface and detect sharp features. Returns maximum
// curvature (per surface point. Note: returns per point, not
// per localpoint)
const auto& ts = *tsPtr;
auto tcurv(triSurfaceTools::curvatures(ts));
auto& curv = tcurv.ref();
// Reset curvature on sharp edges (and neighbours since
// curvature uses extended stencil)
{
const auto& edgeFaces = ts.edgeFaces();
const auto& edges = ts.edges();
const auto& points = ts.points();
const auto& mp = ts.meshPoints();
bitSet isOnSharpEdge(points.size());
forAll(edgeFaces, edgei)
{
const auto& eFaces = edgeFaces[edgei];
const edge meshE(mp, edges[edgei]);
if (eFaces.size() == 2)
{
const auto& f0 = ts[eFaces[0]];
const auto& f1 = ts[eFaces[1]];
const vector n0 = f0.unitNormal(points);
const int dir0 = f0.edgeDirection(meshE);
const int dir1 = f1.edgeDirection(meshE);
vector n1 = f1.unitNormal(points);
if (dir0 == dir1)
{
// Flip since use edge in same direction
// (should not be the case for 'proper'
// surfaces)
n1 = -n1;
}
if ((n0&n1) < cosAngle)
{
isOnSharpEdge.set(meshE[0]);
isOnSharpEdge.set(meshE[1]);
}
}
}
// Extend by one layer
{
bitSet oldOnSharpEdge(isOnSharpEdge);
isOnSharpEdge = false;
for (const auto& f : ts)
{
for (const label pointi : f)
{
if (oldOnSharpEdge[pointi])
{
// Mark all points on triangle
isOnSharpEdge.set(f);
break;
}
}
}
}
// Unmark curvature
autoPtr<OBJstream> str;
//if (debug)
//{
// str.reset
// (
// new OBJstream
// (
// "sharpEdgePoints_"
// + geom.name()
// + ".obj"
// )
// );
// Info<< "Writing sharp edge points to "
// << str().name() << endl;
//}
for (const label pointi : isOnSharpEdge)
{
// Reset measured point-based curvature
curv[pointi] = 0.0;
if (str)
{
str().write(points[pointi]);
}
}
}
// Reset curvature on -almost- sharp edges.
// This resets the point-based curvature if the edge
// is considered to be a sharp edge based on its actual
// curvature. This is only used if the 'ignore' level is
// given.
{
// Pass 1: accumulate constraints on the points - get
// the minimum of curvature constraints on the
// connected triangles. Looks at user-specified
// min curvature - does not use actual measured
// curvature
scalarField pointMinCurv(ts.nPoints(), VGREAT);
forAll(ts, i)
{
// Is ignore level given for surface
const label level = curvIgnore[i];
if (level >= 0)
{
// Convert to (inv) size
const scalar length = level0EdgeLength/(2<<level);
const scalar invLength = 1.0/length;
for (const label pointi : ts[i])
{
if
(
invLength < pointMinCurv[pointi]
&& curv[pointi] > SMALL
)
{
//Pout<< "** at location:"
// << ts.points()[pointi]
// << " measured curv:" << curv[pointi]
// << " radius:" << 1.0/curv[pointi]
// << " ignore level:" << level
// << " ignore radius:" << length
// << " resetting minCurv to "
// << invLength
// << endl;
}
pointMinCurv[pointi] =
min(pointMinCurv[pointi], invLength);
}
}
}
// Clip curvature (should do nothing for most points unless
// ignore-level is triggered)
forAll(pointMinCurv, pointi)
{
if (pointMinCurv[pointi] < curv[pointi])
{
//Pout<< "** at location:" << ts.points()[pointi]
// << " measured curv:" << curv[pointi]
// << " radius:" << 1.0/curv[pointi]
// << " cellLimit:" << pointMinCurv[pointi]
// << endl;
// Set up to ignore point
//curv[pointi] = pointMinCurv[pointi];
curv[pointi] = 0.0;
}
}
}
forAll(ts, i)
{
const auto& f = ts[i];
// Take max curvature (= min radius of curvature)
cellCurv[i] = max(curv[f[0]], max(curv[f[1]], curv[f[2]]));
}
}
//if(debug)
//{
// const scalar maxCurv = gMax(cellCurv);
// if (maxCurv > SMALL)
// {
// const scalar r = scalar(1.0)/maxCurv;
//
// Pout<< "For geometry " << geom.name()
// << " have curvature max " << maxCurv
// << " which equates to radius:" << r
// << " which equates to refinement level "
// << log2(level0EdgeLength/r)
// << endl;
// }
//}
forAll(cellCurv, i)
{
if (cellCurv[i] > SMALL && nCurvCells[i] > 0)
{
//- ?If locally have a cached field override the
// surface-based level ignore any curvature?
//if (minLevelField[i] > surfMin[i])
//{
// // Ignore curvature
//}
//else
//if (surfMin[i] == surfMax[i])
//{
// // Ignore curvature. Bypass calculation below.
//}
//else
{
// Re-work the curvature into a radius and into a
// number of cells
const scalar r = scalar(1.0)/cellCurv[i];
const scalar level =
log2(nCurvCells[i]*level0EdgeLength/r);
const label l = round(level);
if (l > minLevelField[i] && l <= surfMax[i])
{
minLevelField[i] = l;
}
}
}
}
// Store minLevelField on surface
const_cast<searchableSurface&>(geom).setField(minLevelField);
}
}
}
// Find intersections of edge. Return -1 or first surface with higher minLevel
// number.
void Foam::refinementSurfaces::findHigherIntersection

29
src/mesh/snappyHexMesh/refinementSurfaces/refinementSurfaces.H
View File

@ -110,6 +110,9 @@ class refinementSurfaces
// (in gap refinement)
boolList selfProximity_;
//- From global region number to curvature specification
List<FixedList<label, 4>> extendedCurvatureLevel_;
//- From global region number to perpendicular angle
scalarField perpendicularAngle_;
@ -263,6 +266,19 @@ public:
return selfProximity_;
}
//- From global region number to specification of curvature
// refinement: 4 labels specifying
// - minimum wanted number of cells in the curvature radius
// - ?minimum cell level when to start trying to detect gaps
// - maximum cell level to refine to (so do not detect curvature if
// celllevel >= maximum level)
// - minimum radius to ignore (expressed as refinement level).
// This can be used to ignore feature-edges. Set to -1 to ignore.
const List<FixedList<label, 4>>& extendedCurvatureLevel() const
{
return extendedCurvatureLevel_;
}
//- From global region number to perpendicular angle
const scalarField& perpendicularAngle() const
{
@ -307,12 +323,23 @@ public:
//- Per surface the maximum extendedGapLevel over all its regions
labelList maxGapLevel() const;
//- Calculate minLevelFields
//- Per surface the maximum curvatureLevel over all its regions
labelList maxCurvatureLevel() const;
//- Calculate minLevelFields according to both surface- and
// shell-based levels
void setMinLevelFields
(
const shellSurfaces& shells
);
//- Update minLevelFields according to (triSurface-only) curvature
void setCurvatureMinLevelFields
(
const scalar cosAngle,
const scalar level0EdgeLength
);
////- Helper: count number of triangles per region
//static labelList countRegions(const triSurface&);

13
src/mesh/snappyHexMesh/snappyHexMeshDriver/snappyRefineDriver.C
View File

@ -230,10 +230,16 @@ Foam::label Foam::snappyRefineDriver::smallFeatureRefine
label iter = 0;
// See if any surface has an extendedGapLevel
labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
const labelList surfaceMaxLevel(meshRefiner_.surfaces().maxGapLevel());
const labelList shellMaxLevel(meshRefiner_.shells().maxGapLevel());
const labelList curvMaxLevel(meshRefiner_.surfaces().maxCurvatureLevel());
if (max(surfaceMaxLevel) == 0 && max(shellMaxLevel) == 0)
if
(
max(surfaceMaxLevel) == 0
&& max(shellMaxLevel) == 0
&& max(curvMaxLevel) == 0
)
{
return iter;
}
@ -3382,6 +3388,7 @@ void Foam::snappyRefineDriver::doRefine
(
max(meshRefiner_.surfaces().maxGapLevel()) > 0
|| max(meshRefiner_.shells().maxGapLevel()) > 0
|| max(meshRefiner_.surfaces().maxCurvatureLevel()) > 0
)
{
// In case we use automatic gap level refinement do some pre-refinement