Problem with `model_to_mesh` in `dolfinx.io.gmshio`

mesh
1

Hello,

I am trying to set a mesh using gmsh, and I am struggling with the interface in dolfinx.io.gmshio. My code is the following :

import gmsh
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI

mesh_comm = MPI.COMM_WORLD
gdim = 2
model_rank = 0
gmsh.initialize()

L = 5.0 
H = 1.0 
hmin = 0.1; # min element size
hmax = 0.01; # max element size

p0 = gmsh.model.occ.addPoint(0, 0, 0, hmin, 1)
p1 = gmsh.model.occ.addPoint(0, L, 0, hmin, 2)
p2 = gmsh.model.occ.addPoint(L, H, 0, hmax, 3)
p3 = gmsh.model.occ.addPoint(0, H, 0, hmax, 4)
lines = []
lines.append(gmsh.model.occ.addLine(1, 2, 1))
lines.append(gmsh.model.occ.addLine(2, 3, 2))
lines.append(gmsh.model.occ.addLine(3, 4, 3))
lines.append(gmsh.model.occ.addLine(4, 1, 4))
loop = gmsh.model.occ.addCurveLoop(lines ,5)
surface = gmsh.model.occ.addPlaneSurface([5],6)
gmsh.model.addPhysicalGroup(2, [surface], 1)
gmsh.model.setPhysicalName(2, 1, "surface")
gmsh.model.occ.synchronize()

gmsh.model.mesh.generate(gdim)
domain, cell_tags, facet_tags = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

gmsh.finalize()

I get the following error:

File ~/anaconda3/envs/fenicsx-0.5.1/lib/python3.9/site-packages/dolfinx/io/gmshio.py:205, in model_to_mesh(model, comm, rank, gdim)
    202 perm_sort = np.argsort(cell_dimensions)
    204 # Broadcast cell type data and geometric dimension
--> 205 cell_id = cell_information[perm_sort[-1]]["id"]
    206 tdim = cell_information[perm_sort[-1]]["dim"]
    207 num_nodes = cell_information[perm_sort[-1]]["num_nodes"]

IndexError: index -1 is out of bounds for axis 0 with size 0

It seems, from a similar topic on this forum, that this has something to do with makers. But I am not sure how to solve this issue.

I used to export the gmsh files in msh format and to convert them into xml using meshconvert.convert2xml() from dolfin_utils.meshconvert, but I’m not sure a similar option exists in dolfinx, so I am stuck there…

Many thanks in advance for your help with this!

Claire

2

I have tried adding some tags, but now the kernel dies when I run the code…

import gmsh
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI

mesh_comm = MPI.COMM_WORLD
gdim = 2
model_rank = 0

gmsh.initialize()

L = 5.0 #5
H = 1.0 #1
hmin = 0.1; # min element size
hmax = 0.01; # max element size


p0 = gmsh.model.occ.addPoint(0, 0, 0, hmin, 1)
p1 = gmsh.model.occ.addPoint(0, L, 0, hmin, 2)
p2 = gmsh.model.occ.addPoint(L, H, 0, hmax, 3)
p3 = gmsh.model.occ.addPoint(0, H, 0, hmax, 4)
lines = []
lines.append(gmsh.model.occ.addLine(1, 2, 1))
lines.append(gmsh.model.occ.addLine(2, 3, 2))
lines.append(gmsh.model.occ.addLine(3, 4, 3))
lines.append(gmsh.model.occ.addLine(4, 1, 4))
loop = gmsh.model.occ.addCurveLoop(lines ,5)
surface = gmsh.model.occ.addPlaneSurface([5],6)
gmsh.model.occ.synchronize()

# Add physical tag for bulk
gmsh.model.addPhysicalGroup(2, [surface], 1)
gmsh.model.setPhysicalName(2, 1, "surface")

# Add physical tag for boundaries
gmsh.model.addPhysicalGroup(1, [lines[0]], 1)
gmsh.model.addPhysicalGroup(1, [lines[1]], 2)
gmsh.model.addPhysicalGroup(1, [lines[2]], 3)
gmsh.model.addPhysicalGroup(1, [lines[3]], 4)

gmsh.model.mesh.generate(gdim)
domain, cell_tags, facet_tags = model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)

gmsh.finalize()
3

I guess I’m a buit to late to this post. The issue is the following lines:

This happens because GMSH doesn’t use these lines in the model when one create curve loops, and thus you end up with duplicate intervals.

Here is a minimal example:

from dolfinx.io import XDMFFile
import gmsh
from dolfinx.io.gmshio import model_to_mesh
from mpi4py import MPI

mesh_comm = MPI.COMM_WORLD
gdim = 2
model_rank = 0

gmsh.initialize()

L = 5.0  # 5
H = 1.0  # 1
hmin = 0.1  # min element size
hmax = 0.01  # max element size


p0 = gmsh.model.occ.addPoint(0, 0, 0, hmin, 1)
p1 = gmsh.model.occ.addPoint(0, L, 0, hmin, 2)
p2 = gmsh.model.occ.addPoint(L, H, 0, hmax, 3)
p3 = gmsh.model.occ.addPoint(0, H, 0, hmax, 4)
lines = []
lines.append(gmsh.model.occ.addLine(1, 2, 1))
lines.append(gmsh.model.occ.addLine(2, 3, 2))
lines.append(gmsh.model.occ.addLine(3, 4, 3))
lines.append(gmsh.model.occ.addLine(4, 1, 4))
loop = gmsh.model.occ.addCurveLoop(lines, 5)
surface = gmsh.model.occ.addPlaneSurface([5], 6)
gmsh.model.occ.synchronize()

# Add physical tag for bulk
gmsh.model.addPhysicalGroup(2, [surface], 1)
gmsh.model.setPhysicalName(2, 1, "surface")

# Add physical tag for boundaries
boundary = gmsh.model.getBoundary(
    [(2, surface)], combined=False, oriented=False, recursive=False)
for i, (boundary) in enumerate(boundary):
    gmsh.model.addPhysicalGroup(boundary[0], [boundary[1]], i)

gmsh.model.mesh.generate(gdim)
domain, cell_tags, facet_tags = model_to_mesh(
    gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.finalize()

with XDMFFile(domain.comm, "domain.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_meshtags(cell_tags, domain.geometry)
    xdmf.write_meshtags(facet_tags, domain.geometry)