Designing a simple mesh with boundary layers
I am working to develop a cylindrical mesh with refinement specifically along the longitudinal axis, i.e. with boundary layers in the z direction to mimic meshes like this generated in COMSOL.
Currently, I am using the gmsh Python API, but am having some issues getting the exact mesh I hope for.
I didn't have success using "extrude" or "extrude_boundary_layer", so I ended up using "fragment" in the occ kernel to embed disks in the cylinder, forcing the mesh to conform with these structures. However, the resulting boundary layers do not extend to the edges of the geometry, as shown. (this is because the embedded 2D objects cannot intersect with the boundary of the cylinder when using "fragment")
I appreciate any input on how to better implement this mesh; it may be easily accomplished with "extrudeBoundaryLayer", so please let me know if you can provide guidance on use of that function. My current function is included below:
import pathlib
import gmsh
from typing import Tuple
def create_cylinder_zLayers(
radius: float = 1.0,
height: float = 10.0,
boundary_layers: Tuple = (0.0, 0.0),
h_fine: float = 0.0,
h_coarse: float = 0.0,
surf_tag: int = 10,
vol_tag: int = 1):
import gmsh
# assign refinement parameters if not set
if h_fine == 0.0:
h_fine = 0.08 * radius
if h_coarse == 0.0:
h_coarse = 0.6 * radius
if boundary_layers[0] == 0.0:
boundary_layers = (0.005 * radius, 1.5)
# Create the cylinder using gmsh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)
gmsh.model.add("my_cylinder")
cyl_tag = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, height, radius)
r = boundary_layers[1] # each boundary layer is thicker than the last by this factor
z_vals = [boundary_layers[0]] # thickness of first layer
i = 1
cur_thickness = z_vals[0] * r**i
# create a list of disks to embed in the cylinder to enforce boundary layer refinement
# terminate this process when the thickness is greater than h_coarse or z_cur < 0
while cur_thickness < h_coarse:
z_vals.append(z_vals[-1] + cur_thickness)
z_cur = height - z_vals[-1]
if z_cur <= 0:
break
else:
disk_cur = gmsh.model.occ.add_disk(0, 0, z_cur, 0.95*radius, 0.95*radius)
my_cyl = gmsh.model.occ.fragment([(3, cyl_tag)], [(2, disk_cur)])
cyl_tag = my_cyl[0][0][1]
# update thickness
i = i + 1
cur_thickness = z_vals[0] * r**i
gmsh.model.occ.synchronize()
# assign physical groups - top of cylinder is marked with surf_tag
# and the cylinder volume is marked with vol_tag
gmsh.model.add_physical_group(3, [cyl_tag], tag=vol_tag)
facets = gmsh.model.getBoundary([(3, cyl_tag)])
gmsh.model.add_physical_group(2, [facets[1][1]], tag=surf_tag)
def meshSizeCallback(dim, tag, x, y, z, lc):
# return smaller mesh size when z > 0.8*height
thresh = 0.8
if z > height * thresh:
return h_fine + (h_coarse - h_fine) * (height - z) / (height * (1 - thresh))
else:
return h_coarse
gmsh.model.mesh.setSizeCallback(meshSizeCallback)
# set off the other options for mesh size determination
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 1)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
# this changes the algorithm from Frontal-Delaunay to Delaunay,
# which may provide better results when there are larger gradients in mesh size
gmsh.option.setNumber("Mesh.Algorithm", 5)
gmsh.model.mesh.generate(3)
tmp_folder = pathlib.Path(f"cylinder_test")
tmp_folder.mkdir(exist_ok=True)
gmsh_file = tmp_folder / "cylinder.msh"
gmsh.write(str(gmsh_file)) # save locally
gmsh.finalize()
create_cylinder_zLayers(radius = 0.2, height = 0.6)