I am constructing a Backward Facing Step mesh using gmsh(4.12.1).
This is my script to create the .msh file:
import gmsh
import numpy as np
from mpi4py import MPI
gmsh.initialize()
proc = MPI.COMM_WORLD.rank
if proc == 0:
gmsh.model.add("two_jets")
# Parameters
inf_width = 0.2
step_height = 0.1
step_dist = 1.0
top_len = 3.0
control_width = 0.01
# Mesh sizes
inflow_size = 0.02
control_size = 0.004
step_size = 0.016
outflow_size = 0.03
# Derived parameters
outflow_width = step_height + inf_width
# Points
points = [
(0, step_height, 0, inflow_size), # p1
(step_dist, step_height, 0, step_size), # p2
(step_dist, (step_height * 3/4) + control_width, 0, control_size), # p3
(step_dist, step_height * 3/4, 0, control_size), # p4
(step_dist, step_height * 1/4, 0, control_size), # p5
(step_dist, (step_height * 1/4) - control_width, 0, control_size), # p6
(step_dist, 0, 0, step_size), # p7
(top_len, 0, 0, outflow_size), # p8
(top_len, outflow_width, 0, outflow_size), # p9
(0, outflow_width, 0, inflow_size) # p10
]
ptags = []
for i, (x, y, z, s) in enumerate(points):
ptags.append(gmsh.model.geo.addPoint(x, y, z, s))
# Lines
lines = [
(ptags[0], ptags[1]), #l0
(ptags[1], ptags[2]), #l1
(ptags[2], ptags[3]), #l2
(ptags[3], ptags[4]), #l3
(ptags[4], ptags[5]), #l4
(ptags[5], ptags[6]), #l5
(ptags[6], ptags[7]), #l6
(ptags[7], ptags[8]), #l7
(ptags[8], ptags[9]), #l8
(ptags[9], ptags[0]) #l9
]
ltags = []
for start, end in lines:
ltags.append(gmsh.model.geo.addLine(start, end))
# Curve loop and surface
loop = gmsh.model.geo.addCurveLoop(ltags)
surface = gmsh.model.geo.addPlaneSurface([loop])
gmsh.model.geo.synchronize()
# Physical Groups
gmsh.model.addPhysicalGroup(1, [ltags[9]], 1) # Inlet
gmsh.model.setPhysicalName(1, 1, "Inflow")
gmsh.model.addPhysicalGroup(1, [ltags[7]], 2) # Outlet
gmsh.model.setPhysicalName(1, 2, "Outflow")
gmsh.model.addPhysicalGroup(1, [ltags[0]], 3) # Wall1
gmsh.model.setPhysicalName(1, 3, "Wall1")
gmsh.model.addPhysicalGroup(1, [ltags[1]], 4) # Wall2
gmsh.model.setPhysicalName(1, 4, "Wall2")
gmsh.model.addPhysicalGroup(1, [ltags[3]], 5) # Wall3
gmsh.model.setPhysicalName(1, 5, "Wall3")
gmsh.model.addPhysicalGroup(1, [ltags[5]], 8) # Wall4
gmsh.model.setPhysicalName(1, 8, "Wall4")
gmsh.model.addPhysicalGroup(1, [ltags[6]], 9) # Bottom wall
gmsh.model.setPhysicalName(1, 9, "Wall5")
gmsh.model.addPhysicalGroup(1, [ltags[8]], 10) # Bottom wall
gmsh.model.setPhysicalName(1, 10, "Wall6")
gmsh.model.addPhysicalGroup(1, [ltags[2]], 6) # Jet1
gmsh.model.setPhysicalName(1, 6, "Jet1")
gmsh.model.addPhysicalGroup(1, [ltags[4]], 7) # Jet2
gmsh.model.setPhysicalName(1, 7, "Jet2")
gmsh.model.addPhysicalGroup(2, [surface], 5) # Fluid domain
gmsh.model.setPhysicalName(2, 5, "Fluid")
gmsh.model.mesh.generate(2)
gmsh.write("two_jets.msh")
for tag in ltags:
num = len(gmsh.model.getEntitiesForPhysicalGroup(1, tag))
print(f"Curve tag {tag} → mesh entities: {num}")
gmsh.finalize()
then i go to .h5 file with this script using meshio(5.3.5):
import meshio
import numpy as np
from dolfin import *
from mpi4py import MPI
proc = MPI.COMM_WORLD.rank
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={"facet_tags": [cell_data.astype(np.int32)]}
)
return out_mesh
if proc == 0:
msh = meshio.read("two_jets.msh")
# Create clearly named files for mesh (cells) and facets (lines)
triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
line_mesh = create_mesh(msh, "line", prune_z=True)
meshio.write("two_jets_mesh.xdmf", triangle_mesh)
meshio.write("two_jets_facets.xdmf", line_mesh)
MPI.COMM_WORLD.barrier()
mesh = Mesh()
with XDMFFile(MPI.COMM_WORLD, "two_jets_mesh.xdmf") as xdmf:
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1)
with XDMFFile(MPI.COMM_WORLD, "two_jets_facets.xdmf") as xdmf:
xdmf.read(mvc, "facet_tags")
surfaces = MeshFunction("size_t", mesh, mvc)
h5 = HDF5File(MPI.COMM_WORLD, "wall_mesh.h5", "w")
h5.write(mesh, "mesh")
h5.write(surfaces, "facet")
print("Tags:", set(surfaces.array()))
when i do: print("Tags:", set(surfaces.array())) i get Tags: {1, 2, 3, 4, 5, 6, 18446744073709551615, 7, 8, 9, 10}
but if i inspect the physical groups of the .msh file everything seems ok: Line physical tags: [ 1 2 3 4 5 6 7 8 9 10]
where i get the tags doing:
import meshio
import numpy as np
# Load mesh
msh = meshio.read("two_jets.msh")
points = msh.points[:, :2]
line_cells = msh.get_cells_type("line")
line_data = msh.get_cell_data("gmsh:physical", "line")
print("Line data dtype:", line_data.dtype)
print("Line physical tags:", np.unique(line_data))
# Find untagged line facets
untagged_mask = line_data == -1
untagged_lines = line_cells[untagged_mask]
# Plot untagged lines in red
for line in untagged_lines:
coords = points[line]
print(coords)
At first i thought that the problem was with an untagged portion of the boundary, but i cannot find it.
I cannot find the error in the conversion routine, can someone help me?
Thank You!!