Skip to content

GitLab

Commit 6da9df71 authored by mattijs's avatar mattijs
Browse files

ENH: snappyHexMesh: allow layer thickness specification by total thickness

parent 389170bc
Hide whitespace changes
Inline Side-by-side
applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
View file @ 6da9df71
......@@ -226,37 +226,37 @@ castellatedMeshControls
// Settings for the snapping.
snapControls
{
//- Number of patch smoothing iterations before finding correspondence
// to surface
// Number of patch smoothing iterations before finding correspondence
// to surface
nSmoothPatch 3;
//- Maximum relative distance for points to be attracted by surface.
// True distance is this factor times local maximum edge length.
// Maximum relative distance for points to be attracted by surface.
// True distance is this factor times local maximum edge length.
// Note: changed(corrected) w.r.t 17x! (17x used 2* tolerance)
tolerance 2.0;
//- Number of mesh displacement relaxation iterations.
// Number of mesh displacement relaxation iterations.
nSolveIter 30;
//- Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
// Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
nRelaxIter 5;
// Feature snapping
//- Number of feature edge snapping iterations.
// Leave out altogether to disable.
// Number of feature edge snapping iterations.
// Leave out altogether to disable.
nFeatureSnapIter 10;
//- Detect (geometric only) features by sampling the surface
// (default=false).
// Detect (geometric only) features by sampling the surface
// (default=false).
implicitFeatureSnap false;
//- Use castellatedMeshControls::features (default = true)
// Use castellatedMeshControls::features (default = true)
explicitFeatureSnap true;
//- Detect features between multiple surfaces
// (only for explicitFeatureSnap, default = false)
// Detect features between multiple surfaces
// (only for explicitFeatureSnap, default = false)
multiRegionFeatureSnap false;
}
......@@ -267,9 +267,43 @@ addLayersControls
// size of the refined cell outside layer (true) or absolute sizes (false).
relativeSizes true;
// Layer thickness specification. This can be specified in one of four ways
// - expansionRatio and finalLayerThickness (cell nearest internal mesh)
// - expansionRatio and firstLayerThickness (cell on surface)
// - overall thickness and firstLayerThickness
// - overall thickness and finalLayerThickness
// Expansion factor for layer mesh
expansionRatio 1.0;
// Wanted thickness of the layer furthest away from the wall.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
finalLayerThickness 0.3;
// Wanted thickness of the layer next to the wall.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
//firstLayerThickness 0.3;
// Wanted overall thickness of layers.
// If relativeSizes this is relative to undistorted size of cell
// outside layer.
//thickness 0.5
// Minimum overall thickness of total layers. If for any reason layer
// cannot be above minThickness do not add layer.
// If relativeSizes this is relative to undistorted size of cell
// outside layer..
minThickness 0.25;
// Per final patch (so not geometry!) the layer information
// Note: This behaviour changed after 21x. Any non-mentioned patches
// now slide unless nSurfaceLayers is explicitly mentioned to be 0.
// now slide unless:
// - nSurfaceLayers is explicitly mentioned to be 0.
// - angle to nearest surface < slipFeatureAngle (see below)
layers
{
sphere.stl_firstSolid
......@@ -281,9 +315,9 @@ addLayersControls
{
nSurfaceLayers 1;
// Per patch layer data
expansionRatio 1.3;
expansionRatio 1.3;
finalLayerThickness 0.3;
minThickness 0.1;
minThickness 0.1;
}
// Disable any mesh shrinking and layer addition on any point of
......@@ -294,41 +328,24 @@ addLayersControls
}
}
// Expansion factor for layer mesh
expansionRatio 1.0;
//- Wanted thickness of final added cell layer. If multiple layers
// is the
// thickness of the layer furthest away from the wall.
// Relative to undistorted size of cell outside layer.
// is the thickness of the layer furthest away from the wall.
// See relativeSizes parameter.
finalLayerThickness 0.3;
//- Minimum thickness of cell layer. If for any reason layer
// cannot be above minThickness do not add layer.
// Relative to undistorted size of cell outside layer.
// See relativeSizes parameter.
minThickness 0.25;
//- If points get not extruded do nGrow layers of connected faces that are
// also not grown. This helps convergence of the layer addition process
// close to features.
// If points get not extruded do nGrow layers of connected faces that are
// also not grown. This helps convergence of the layer addition process
// close to features.
// Note: changed(corrected) w.r.t 17x! (didn't do anything in 17x)
nGrow 0;
// Advanced settings
//- When not to extrude surface. 0 is flat surface, 90 is when two faces
// make straight angle.
// When not to extrude surface. 0 is flat surface, 90 is when two faces
// make straight angle.
featureAngle 60;
//- At non-patched sides allow mesh to slip if extrusion direction makes
// angle larger than slipFeatureAngle.
// At non-patched sides allow mesh to slip if extrusion direction makes
// angle larger than slipFeatureAngle.
slipFeatureAngle 30;
//- Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
// Maximum number of snapping relaxation iterations. Should stop
// before upon reaching a correct mesh.
nRelaxIter 5;
// Number of smoothing iterations of surface normals
......@@ -374,51 +391,51 @@ addLayersControls
// where to undo.
meshQualityControls
{
//- Maximum non-orthogonality allowed. Set to 180 to disable.
// Maximum non-orthogonality allowed. Set to 180 to disable.
maxNonOrtho 65;
//- Max skewness allowed. Set to <0 to disable.
// Max skewness allowed. Set to <0 to disable.
maxBoundarySkewness 20;
maxInternalSkewness 4;
//- Max concaveness allowed. Is angle (in degrees) below which concavity
// is allowed. 0 is straight face, <0 would be convex face.
// Set to 180 to disable.
// Max concaveness allowed. Is angle (in degrees) below which concavity
// is allowed. 0 is straight face, <0 would be convex face.
// Set to 180 to disable.
maxConcave 80;
//- Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
// Minimum pyramid volume. Is absolute volume of cell pyramid.
// Set to a sensible fraction of the smallest cell volume expected.
// Set to very negative number (e.g. -1E30) to disable.
minVol 1e-13;
//- Minimum quality of the tet formed by the face-centre
// and variable base point minimum decomposition triangles and
// the cell centre. This has to be a positive number for tracking
// to work. Set to very negative number (e.g. -1E30) to
// disable.
// <0 = inside out tet,
// 0 = flat tet
// 1 = regular tet
// Minimum quality of the tet formed by the face-centre
// and variable base point minimum decomposition triangles and
// the cell centre. This has to be a positive number for tracking
// to work. Set to very negative number (e.g. -1E30) to
// disable.
// <0 = inside out tet,
// 0 = flat tet
// 1 = regular tet
minTetQuality 1e-9;
//- Minimum face area. Set to <0 to disable.
// Minimum face area. Set to <0 to disable.
minArea -1;
//- Minimum face twist. Set to <-1 to disable. dot product of face normal
//- and face centre triangles normal
// Minimum face twist. Set to <-1 to disable. dot product of face normal
// and face centre triangles normal
minTwist 0.05;
//- minimum normalised cell determinant
//- 1 = hex, <= 0 = folded or flattened illegal cell
// minimum normalised cell determinant
// 1 = hex, <= 0 = folded or flattened illegal cell
minDeterminant 0.001;
//- minFaceWeight (0 -> 0.5)
// minFaceWeight (0 -> 0.5)
minFaceWeight 0.05;
//- minVolRatio (0 -> 1)
// minVolRatio (0 -> 1)
minVolRatio 0.01;
//must be >0 for Fluent compatibility
// must be >0 for Fluent compatibility
minTriangleTwist -1;
//- if >0 : preserve single cells with all points on the surface if the
......@@ -429,9 +446,9 @@ meshQualityControls
// Advanced
//- Number of error distribution iterations
// Number of error distribution iterations
nSmoothScale 4;
//- amount to scale back displacement at error points
// amount to scale back displacement at error points
errorReduction 0.75;
// Optional : some meshing phases allow usage of relaxed rules.
......
src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
View file @ 6da9df71
......@@ -318,7 +318,6 @@ void Foam::autoLayerDriver::handleNonManifolds
}
}
Info<< "Set displacement to zero for all " << nNonManif
<< " non-manifold points" << endl;
}
......@@ -443,86 +442,6 @@ void Foam::autoLayerDriver::handleFeatureAngle
}
//Foam::tmp<Foam::scalarField> Foam::autoLayerDriver::undistortedEdgeLength
//(
// const indirectPrimitivePatch& pp,
// const bool relativeSizes,
// const bool faceSize
//)
//{
// const fvMesh& mesh = meshRefiner_.mesh();
//
// tmp<scalarField> tfld(new scalarField());
// scalarField& fld = tfld();
//
// if (faceSize)
// {
// fld.setSize(pp.size());
// }
// else
// {
// fld.setSize(pp.nPoints());
// }
//
//
// if (relativeSizes)
// {
// const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength();
// const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel();
//
// if (faceSize)
// {
// forAll(pp, i)
// {
// label faceI = pp.addressing()[i];
// label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
// fld[i] = edge0Len/(1<<ownLevel);
// }
// }
// else
// {
// // Determine per point the max cell level of connected cells
// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//
// labelList maxPointLevel(pp.nPoints(), labelMin);
//
// forAll(pp, i)
// {
// label faceI = pp.addressing()[i];
// label ownLevel = cellLevel[mesh.faceOwner()[faceI]];
//
// const face& f = pp.localFaces()[i];
// forAll(f, fp)
// {
// maxPointLevel[f[fp]] =
// max(maxPointLevel[f[fp]], ownLevel);
// }
// }
//
// syncTools::syncPointList
// (
// mesh,
// pp.meshPoints(),
// maxPointLevel,
// maxEqOp<label>(),
// labelMin // null value
// );
//
//
// forAll(maxPointLevel, pointI)
// {
// // Find undistorted edge size for this level.
// fld[i] = edge0Len/(1<<maxPointLevel[pointI]);
// }
// }
// }
// else
// {
// // Use actual cell size
// }
//}
// No extrusion on cells with warped faces. Calculates the thickness of the
// layer and compares it to the space the warped face takes up. Disables
// extrusion if layer thickness is more than faceRatio of the thickness of
......@@ -702,7 +621,6 @@ void Foam::autoLayerDriver::handleWarpedFaces
//}
// No extrusion on faces with differing number of layers for points
void Foam::autoLayerDriver::setNumLayers
(
const labelList& patchToNLayers,
......@@ -865,15 +783,6 @@ Foam::autoLayerDriver::makeLayerDisplacementField
}
//Pout<< "*** makeLayerDisplacementField : boundary conditions:" << endl;
//forAll(patchFieldTypes, patchI)
//{
// Pout<< "\t" << patchI << " name:" << pointPatches[patchI].name()
// << " type:" << patchFieldTypes[patchI]
// << " nLayers:" << numLayers[patchI]
// << endl;
//}
const polyMesh& mesh = pMesh();
// Note: time().timeName() instead of meshRefinement::timeName() since
......@@ -1053,11 +962,10 @@ void Foam::autoLayerDriver::determineSidePatches
patchDict.add("nFaces", 0);
patchDict.add("startFace", mesh.nFaces());
Pout<< "Adding patch " << patchI
<< " name:" << name
<< " between " << Pstream::myProcNo()
<< " and " << nbrProcI
<< endl;
//Pout<< "Adding patch " << patchI
// << " name:" << name
// << " between " << Pstream::myProcNo()
// << " and " << nbrProcI << endl;
label procPatchI = meshRefiner_.appendPatch
(
......@@ -1090,12 +998,7 @@ void Foam::autoLayerDriver::calculateLayerThickness
(
const indirectPrimitivePatch& pp,
const labelList& patchIDs,
const scalarField& patchExpansionRatio,
const bool relativeSizes,
const scalarField& patchFinalLayerThickness,
const scalarField& patchMinThickness,
const layerParameters& layerParams,
const labelList& cellLevel,
const labelList& patchNLayers,
const scalar edge0Len,
......@@ -1111,12 +1014,13 @@ void Foam::autoLayerDriver::calculateLayerThickness
// Rework patch-wise layer parameters into minimum per point
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Note: only layer parameters consistent with layer specification
// method (see layerParameters) will be correct.
scalarField firstLayerThickness(pp.nPoints(), GREAT);
scalarField finalLayerThickness(pp.nPoints(), GREAT);
scalarField totalThickness(pp.nPoints(), GREAT);
scalarField expRatio(pp.nPoints(), GREAT);
// Reuse input fields
expansionRatio.setSize(pp.nPoints());
expansionRatio = GREAT;
thickness.setSize(pp.nPoints());
thickness = GREAT;
minThickness.setSize(pp.nPoints());
minThickness = GREAT;
......@@ -1130,20 +1034,30 @@ void Foam::autoLayerDriver::calculateLayerThickness
{
label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
expansionRatio[ppPointI] = min
firstLayerThickness[ppPointI] = min
(
expansionRatio[ppPointI],
patchExpansionRatio[patchI]
firstLayerThickness[ppPointI],
layerParams.firstLayerThickness()[patchI]
);
thickness[ppPointI] = min
finalLayerThickness[ppPointI] = min
(
thickness[ppPointI],
patchFinalLayerThickness[patchI]
finalLayerThickness[ppPointI],
layerParams.finalLayerThickness()[patchI]
);
totalThickness[ppPointI] = min
(
totalThickness[ppPointI],
layerParams.thickness()[patchI]
);
expRatio[ppPointI] = min
(
expRatio[ppPointI],
layerParams.expansionRatio()[patchI]
);
minThickness[ppPointI] = min
(
minThickness[ppPointI],
patchMinThickness[patchI]
layerParams.minThickness()[patchI]
);
}
}
......@@ -1152,7 +1066,7 @@ void Foam::autoLayerDriver::calculateLayerThickness
(
mesh,
pp.meshPoints(),
expansionRatio,
firstLayerThickness,
minEqOp<scalar>(),
GREAT // null value
);
......@@ -1160,7 +1074,23 @@ void Foam::autoLayerDriver::calculateLayerThickness
(
mesh,
pp.meshPoints(),
thickness,
finalLayerThickness,
minEqOp<scalar>(),
GREAT // null value
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
totalThickness,
minEqOp<scalar>(),
GREAT // null value
);
syncTools::syncPointList
(
mesh,
pp.meshPoints(),
expRatio,
minEqOp<scalar>(),
GREAT // null value
);
......@@ -1182,14 +1112,18 @@ void Foam::autoLayerDriver::calculateLayerThickness
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// by multiplying with the internal cell size.
if (relativeSizes)
if (layerParams.relativeSizes())
{
if (min(patchMinThickness) < 0 || max(patchMinThickness) > 2)
if
(
min(layerParams.minThickness()) < 0
|| max(layerParams.minThickness()) > 2
)
{
FatalErrorIn("calculateLayerThickness(..)")
<< "Thickness should be factor of local undistorted cell size."
<< " Valid values are [0..2]." << nl
<< " minThickness:" << patchMinThickness
<< " minThickness:" << layerParams.minThickness()
<< exit(FatalError);
}
......@@ -1225,38 +1159,114 @@ void Foam::autoLayerDriver::calculateLayerThickness
{
// Find undistorted edge size for this level.
scalar edgeLen = edge0Len/(1<<maxPointLevel[pointI]);
thickness[pointI] *= edgeLen;
firstLayerThickness[pointI] *= edgeLen;
finalLayerThickness[pointI] *= edgeLen;
totalThickness[pointI] *= edgeLen;
minThickness[pointI] *= edgeLen;
}
}
// Rework thickness (of final layer) into overall thickness of all layers
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Rework thickness parameters into overall thickness
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
forAll(thickness, pointI)
forAll(firstLayerThickness, pointI)
{
// Calculate layer thickness based on expansion ratio
// and final layer height
if (expansionRatio[pointI] == 1)
thickness[pointI] = layerParams.layerThickness
(
patchNLayers[pointI],
firstLayerThickness[pointI],
finalLayerThickness[pointI],
totalThickness[pointI],
expRatio[pointI]
);
expansionRatio[pointI] = layerParams.layerExpansionRatio
(
patchNLayers[pointI],
firstLayerThickness[pointI],
finalLayerThickness[pointI],
totalThickness[pointI],
expRatio[pointI]
);
}
//Info<< "calculateLayerThickness : min:" << gMin(thickness)
// << " max:" << gMax(thickness) << endl;
// Print a bit
{
const polyBoundaryMesh& patches = mesh.boundaryMesh();
// Find maximum length of a patch name, for a nicer output
label maxPatchNameLen = 0;
forAll(patchIDs, i)
{
thickness[pointI] *= patchNLayers[pointI];
label patchI = patchIDs[i];
word patchName = patches[patchI].name();
maxPatchNameLen = max(maxPatchNameLen, label(patchName.size()));
}
else
Info<< nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << "patch"
<< setw(0) << " faces layers avg thickness[m]" << nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << " "
<< setw(0) << " near-wall overall" << nl
<< setf(ios_base::left) << setw(maxPatchNameLen) << "-----"
<< setw(0) << " ----- ------ --------- -------" << endl;
forAll(patchIDs, i)
{
label patchI = patchIDs[i];
scalar invExpansion = 1.0 / expansionRatio[pointI];
label nLay = patchNLayers[pointI];
thickness[pointI] *=
(1.0 - pow(invExpansion, nLay))
/ (1.0 - invExpansion);
}
}
const labelList& meshPoints = patches[patchI].meshPoints();
scalar sumThickness = 0;
scalar sumNearWallThickness = 0;
//Info<< "calculateLayerThickness : min:" << gMin(thickness)
// << " max:" << gMax(thickness) << endl;
forAll(meshPoints, patchPointI)
{
label ppPointI = pp.meshPointMap()[meshPoints[patchPointI]];
sumThickness += thickness[ppPointI];
sumNearWallThickness += layerParams.firstLayerThickness
(
patchNLayers[ppPointI],
firstLayerThickness[ppPointI],
finalLayerThickness[ppPointI],
thickness[ppPointI],
expansionRatio[ppPointI]
);
}
label totNPoints = returnReduce(meshPoints.size(), sumOp<label>());
// For empty patches, totNPoints is 0.
scalar avgThickness = 0;