Hi all,
I have an issue when generating a boundary layer in GMSH. I am making a mesh with a boundary layer perpendicular to a boundary, as seen here. Overall it is doing a pretty good job.
The issue is at the corner. It seems like the boundary layer extends into the outer domain, even though "ExcludedSurfacesList" includes the outer domain. Here is a zoomed in corner.
I tried asking on the gmsh gitlab, but it requires a gitlab account and I had difficulties signing up as detailed here
So I really hope here someone might know how to help!
Many thanks in advance,
Katie
Here is a mwe
import gmsh
import meshio
from mpi4py import MPI
from dolfinx import io
gmsh.initialize()
lc = 0.1
gmsh.model.geo.addPoint(0.0, -5, 0, lc, 1)
gmsh.model.geo.addPoint(3, -5, 0, lc, 2)
gmsh.model.geo.addPoint(3, 5, 0, lc, 3)
gmsh.model.geo.addPoint(0.0, 5, 0, lc, 4)
gmsh.model.geo.addPoint(0.0, 1, 0, lc, 5)
gmsh.model.geo.addPoint(1, 1, 0, lc, 6)
gmsh.model.geo.addPoint(1, - 1, 0, lc, 7)
gmsh.model.geo.addPoint(0.0, -1, 0, lc, 8)
gmsh.model.geo.addLine(1, 2, 1) # adding curved lines for billet and air domain
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 5, 4)
gmsh.model.geo.addLine(5, 6, 5)
gmsh.model.geo.addLine(6, 7, 6)
gmsh.model.geo.addLine(7, 8, 7)
gmsh.model.geo.addLine(8, 1, 8)
gmsh.model.geo.addLine(5, 8, 9)
gmsh.model.geo.addCurveLoop([1, 2, 3, 4,5,6,7,8], 1)
gmsh.model.geo.addCurveLoop([5,6,7,-9], 2)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.geo.addPlaneSurface([2], 2)
gmsh.model.geo.synchronize()
air_right_boundary = 3
air_topbottom_boundary = 4
air_left_boundary = 5
gmsh.model.addPhysicalGroup(1, [2], air_right_boundary)
gmsh.model.addPhysicalGroup(1, [1,3], air_topbottom_boundary)
gmsh.model.addPhysicalGroup(1, [4,9,8], air_left_boundary)
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.addPhysicalGroup(2, [2], 2)
# Adding field sizes
# Generate mesh
gmsh.model.mesh.field.add("BoundaryLayer", 1)
gmsh.model.mesh.field.setNumbers(1, "CurvesList", [5, 6, 7])
gmsh.model.mesh.field.setNumbers(1, "PointsList", [5,8])
gmsh.model.mesh.field.setNumbers(1, "ExcludedSurfacesList", [1])
gmsh.model.mesh.field.setNumber(1, "SizeFar", 0.1)
gmsh.model.mesh.field.setNumber(1, "Size", 1e-3)
gmsh.model.mesh.field.setNumber(1, "Ratio", 1.1)
gmsh.model.mesh.field.setNumber(1, "Quads", 1)
gmsh.model.mesh.field.setNumber(1, "Thickness", 0.1)
gmsh.model.mesh.field.setNumber(1, "AnisoMax", 0.1)
#gmsh.model.mesh.field.setNumber(1, "IntersectMetrics", 1)
gmsh.model.mesh.field.add("Distance", 3)
gmsh.model.mesh.field.setNumbers(3, "PointsList", [2,3])
gmsh.model.mesh.field.add("Threshold", 4)
gmsh.model.mesh.field.setNumber(4, "InField", 3)
gmsh.model.mesh.field.setNumber(4, "SizeMin", 0.05)
gmsh.model.mesh.field.setNumber(4, "SizeMax", 0.1)
gmsh.model.mesh.field.setNumber(4, "DistMin", 0.2)
gmsh.model.mesh.field.setNumber(4, "DistMax", 1)
gmsh.model.mesh.field.add("Min", 9)
gmsh.model.mesh.field.setNumbers(9, "FieldsList", [1, 4])
gmsh.model.mesh.field.setAsBoundaryLayer(1)
gmsh.model.mesh.field.setAsBackgroundMesh(9)
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.Algorithm", 8) ## Algorithm for gererating quads
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) # subdividing not smoothing
#gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 3) # Algorithm for smoothing quads, 2 or 3
gmsh.model.mesh.generate(2)
# ... and save it to disk
gmsh.write(f"../output/mwe.msh")
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(
points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data]}
)
return out_mesh
msh = meshio.read(f"../output/mwe.msh")
triangle_mesh = create_mesh(msh, "quad", prune_z=True) # prune = True because it is a 2D mesh
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("../output/mesh.xdmf", triangle_mesh) # write the mesh
meshio.write("../output/mt.xdmf", line_mesh)
with io.XDMFFile(MPI.COMM_WORLD, "../output/mesh.xdmf", "r") as xdmf:
whole_domain = xdmf.read_mesh(name="Grid")
ct = xdmf.read_meshtags(whole_domain, name="Grid") # these are the meshtags for the different 2D physical groups!
whole_domain.topology.create_connectivity(whole_domain.topology.dim, whole_domain.topology.dim-1)
whole_domain.topology.create_connectivity(whole_domain.topology.dim-1 , whole_domain.topology.dim)
xdmf_i = io.XDMFFile(whole_domain.comm, "../output/mwe_mesh.xdmf", "w", encoding=io.XDMFFile.Encoding.ASCII)
xdmf_i.write_mesh(whole_domain)
xdmf_i.close()