Error importing mesh with filleted corners

mesh
1

Hi,

I have a simple 2D geometry where I fillet each corner. Everything works in gmsh but somehow it causes the `model_to_mesh` operation to fail. I will appreciate if someone could point a way to get around this. Below is an MWE that was tried on dolfinx 0.9 installed through conda:

import dolfinx, gmsh, numpy as np, mpi4py
from dolfinx.io import gmshio

gmsh.finalize()
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
model_rank = 0

occ = gmsh.model.occ
gdim = 2

fillet = 0.1 # 0.0 means no fillet and the code works
ID = 5
OD = 10
L = 10

pnts_fl = 4*[[]]

pnts_fl[0] = occ.addPoint(ID/2, -L/2, 0)
pnts_fl[1] = occ.addPoint(OD/2, -L/2, 0)
pnts_fl[2] = occ.addPoint(OD/2, L/2, 0)
pnts_fl[3] = occ.addPoint(ID/2, L/2, 0)    

lines_fl = 4*[[]]
for j in range(len(lines_fl)):
    lines_fl[j] = occ.addLine(pnts_fl[j], pnts_fl[(j+1)%len(lines_fl)])

if fillet > 0.0:
    occ.synchronize()
    fillet_tags = 4*[[]]
    for j in range(len(lines_fl)):
        fillet_tags[j] = occ.fillet2D(lines_fl[j], lines_fl[(j+1)%len(lines_fl)], fillet)

    lines_tmp = 8*[[]]
    for j in range(0, 2*len(lines_fl), 2):
        lines_tmp[j] = lines_fl[j//2]
        lines_tmp[j+1] = fillet_tags[j//2]

    lines_fl = lines_tmp
 
loop_fl1 = occ.addCurveLoop(lines_fl)

fl1 = occ.addPlaneSurface([loop_fl1])
occ.synchronize()
gmsh.model.mesh.generate(gdim)

all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

model_rank = 0
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
2

The first check to run is if all nodes loaded in the mesh is part of the cells.

I.e. install meshio, and add the lines:

import gmsh, numpy as np
import meshio

gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0)  # disable output messages
gmsh.clear()
model_rank = 0

occ = gmsh.model.occ
gdim = 2

fillet = 0.1  # 0.0 means no fillet and the code works
ID = 5
OD = 10
L = 10

pnts_fl = 4 * [[]]

pnts_fl[0] = occ.addPoint(ID / 2, -L / 2, 0)
pnts_fl[1] = occ.addPoint(OD / 2, -L / 2, 0)
pnts_fl[2] = occ.addPoint(OD / 2, L / 2, 0)
pnts_fl[3] = occ.addPoint(ID / 2, L / 2, 0)

lines_fl = 4 * [[]]
for j in range(len(lines_fl)):
    lines_fl[j] = occ.addLine(pnts_fl[j], pnts_fl[(j + 1) % len(lines_fl)])

if fillet > 0.0:
    occ.synchronize()
    fillet_tags = 4 * [[]]
    for j in range(len(lines_fl)):
        fillet_tags[j] = occ.fillet2D(
            lines_fl[j], lines_fl[(j + 1) % len(lines_fl)], fillet
        )

    lines_tmp = 8 * [[]]
    for j in range(0, 2 * len(lines_fl), 2):
        lines_tmp[j] = lines_fl[j // 2]
        lines_tmp[j + 1] = fillet_tags[j // 2]

    lines_fl = lines_tmp

loop_fl1 = occ.addCurveLoop(lines_fl)

fl1 = occ.addPlaneSurface([loop_fl1])
occ.synchronize()

all_doms = gmsh.model.getEntities(gdim)
for j, dom in enumerate(all_doms):
    gmsh.model.addPhysicalGroup(dom[0], [dom[1]], j + 1)  # create the main group/node

# number all boundaries
all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(
        edge[0], [edge[1]], edge[1]
    )  # create the main group/node

gmsh.model.mesh.generate(gdim)
gmsh.option.setNumber("Mesh.SaveAll", 0)
gmsh.write("test.msh")
import meshio

mm = meshio.read("test.msh")
lines = mm.cells_dict["line"]
triangles = mm.cells_dict["triangle"]
points = mm.points

# Sanity check 1:
# Check that all points in the mesh is part of the triangular cells
unique_cell_vertices = np.unique(triangles.flatten())
print(points.shape[0], unique_cell_vertices)
assert np.allclose(unique_cell_vertices, np.arange(points.shape[0]))

which returns

90 [ 1  3  5  7  9 11 12 13 16 17 42 43 44 45 46 47 48 49 50 51 52 53 54 55
 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
 80 81 82 83 84 85 86 87 88 89]

i.e. there are 90 nodes read in, but you do not have 90 unique nodes in the cells. This stems from

all_edges = gmsh.model.getEntities(gdim - 1)
for j, edge in enumerate(all_edges):
    gmsh.model.addPhysicalGroup(
        edge[0], [edge[1]], edge[1]
    )  # create the main group/node

where you are not just getting the boundary of your mesh, but also all edges in GMSH, which are not part of your mesh.

Change this to

# number all boundaries
all_edges = gmsh.model.getBoundary(all_doms, oriented=False)
for j, edge in enumerate(all_edges):
   gmsh.model.addPhysicalGroup(edge[0], [edge[1]], edge[1])  # create the main group/node

which then yields the following mesh:

image
image900×1602 152 KB

with


model_rank = 0
from mpi4py import MPI
import dolfinx.io

gmshio = dolfinx.io.gmshio
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_meshtags(ct, mesh.geometry)
    xdmf.write_meshtags(ft, mesh.geometry)
2 Likes
3

I will be able to incorporate your solution tomorrow or during the weekend but for now I cannot thank you enough for such a swift yet detailed reply!