Error extruding gmsh surfaces to abutting volumes

My application deals with many shapes processed with gmsh's Python API, but I have reduced my problem to a MWE of two polygons which lie in separate planes. Their coordinates are:

# polygon spans (x, y, zs[0])
xy1 = [(0.0, 504.0), (247.0, 504.0), (247.0, 496.0), (0.0, 496.0), (0.0, 504.0)]
zs1 = [0, 50]

xy2 = [(0.0, 396.0), (0.0, 496.0), (0.0, 504.0), (247.0, 504.0), (285.0, 504.0), (285.0, 496.0), (285.0, 396.0), (285.0, 388.0), (247.0, 388.0), (0.0, 388.0), (0.0, 396.0)]
zs2 = [50, 100]

My first step is to project these coordinates onto a finite grid so gmsh doesn't complain about node duplication. I do this by merging proximate points, although this isn't necessary for the example above so I simplify it here for clarity:

coord_labels = {}

def get_coord_label(xyz):
    if xyz not in coord_labels:
        coord_labels[xyz] = gmsh.model.geo.add_point(*xyz)
    return coord_labels[xyz]

pts1 = [get_coord_label((x,y,zs1[0])) for (x,y) in xy1]
pts2 = [get_coord_label((x,y,zs2[0])) for (x,y) in xy2]

This produces unique gmsh point labels:

pts1
>>> [1, 2, 3, 4, 1]

pts2
>>> [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 5]

My task now is to extrude these polygons into polyhedra which span zs[0] to zs[1]. Since the polygons are specified in the z = zs[0] plane, I extrude them along vector (0,0,1) by a distance zs[1]-zs[0] like so:

def extrude(pts, dz):
    lines = [gmsh.model.geo.add_line(pts[i], pts[i+1]) for i in range(len(pts)-1)]
    face = gmsh.model.geo.add_plane_surface([gmsh.model.geo.add_curve_loop(lines, reorient=True)])
    objs = gmsh.model.geo.extrude([(2,face)], 0,0, dz)

extrude(pts1, zs1[1] - zs1[0])
extrude(pts2, zs2[1] - zs2[0])

The resulting polyhedra look fine to me when I run

gmsh.model.geo.synchronize()
gmsh.fltk.run()

Above, we see the two volumes; the lower bottom-right one is formed by xy1 and shares its top-surface at z=50 with a sub-bottom-surface of xy2's extrusion.

Alas, this geometry cannot be meshed. Calling

gmsh.model.mesh.generate()

raises exception

No elements in volume 1

If I had to guess, the problem is that gmsh's extrusion function (of xy1) is creating new points within surface z = zs1[1], and are not automatically recognised co-planar with z = zs2[0]. But this is my speculation, and I am not sure why it causes this issue nor how to workaround it.

The entire process (e.g. having to myself group proximate coordinates before extrusion) seems unreasonably difficult for such a simple exercise of extruding polygons! What am I doing wrong? Note in my application, some abutting volumes like this should be merged into a single volume in the mesh, while others should be treated separately (and ergo have an explicit separating surface).

You can use the occ kernel instead of the geo kernel (replace all occurences of "geo" by "occ", and take care with a few function names and arguments). With this change gmsh will mesh both volumes. However the will be meshed seperately because they are not connected in any way. To connect them you need to use the

gmsh.model.occ.fragment(...)

function.

Here is the code:

import gmsh

gmsh.initialize()

# polygon spans (x, y, zs[0])
xy1 = [(0.0, 504.0), (247.0, 504.0), (247.0, 496.0), (0.0, 496.0), (0.0, 504.0)]
zs1 = [0, 50]

xy2 = [(0.0, 396.0), (0.0, 496.0), (0.0, 504.0), (247.0, 504.0), (285.0, 504.0), (285.0, 496.0), (285.0, 396.0), (285.0, 388.0), (247.0, 388.0), (0.0, 388.0), (0.0, 396.0)]
zs2 = [50, 100]

coord_labels = {}

def get_coord_label(xyz):
    if xyz not in coord_labels:
        coord_labels[xyz] = gmsh.model.occ.add_point(*xyz)
    return coord_labels[xyz]


pts1 = [get_coord_label((x,y,zs1[0])) for (x,y) in xy1]
pts2 = [get_coord_label((x,y,zs2[0])) for (x,y) in xy2]
gmsh.model.occ.synchronize()

def extrude(pts, dz):
    lines = [gmsh.model.occ.add_line(pts[i], pts[i+1]) for i in range(len(pts)-1)]
    loop = gmsh.model.occ.addCurveLoop(lines)
    face = gmsh.model.occ.addPlaneSurface([loop])
    objs = gmsh.model.occ.extrude([(2,face)], 0,0, dz)

extrude(pts1, zs1[1] - zs1[0])
extrude(pts2, zs2[1] - zs2[0])
gmsh.model.occ.synchronize()

volTags = gmsh.model.getEntities(3)
gmsh.model.occ.fragment([volTags[0]], [volTags[1]])
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate()
gmsh.fltk.run()
gmsh.finalize()

You should try out commenting the line

gmsh.model.occ.fragment([volTags[0]], [volTags[1]])

and look at the mesh (especially of Volume 2) in the GUI. You will see that the meshes of both volumes do not match at their interface in this case.