Gmsh Python API多边形挤压网格生成报错：No elements in volume 1

我用Gmsh的Python API处理多个几何形状，将问题简化为两个位于不同平面的多边形的最小复现示例，坐标定义如下：

# 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.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]

生成的Gmsh点标签唯一：

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

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

挤压生成多面体 #

接下来我将这些多边形沿(0,0,1)方向挤压成跨越zs[0]到zs[1]的多面体，代码如下：

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])

运行同步和显示代码后，生成的多面体显示正常：

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

其中下方右下的体由xy1生成，其z=50的顶面与xy2挤压体的部分底面重合。但调用网格生成代码：

gmsh.model.mesh.generate()

时抛出异常：

No elements in volume 1

疑问与需求 #

我推测问题在于Gmsh的extrusion函数为xy1生成z=zs1[1]面内的新点，这些点未被自动识别为与z=zs2[0]共面，但不确定原因及解决方法。此外，我认为当前挤压多边形的流程过于繁琐，想了解操作中的错误之处。

注：我的应用场景中，部分相邻体需要合并为单个网格体，其余则需保留分隔面。

内容的提问来源于stack exchange，提问作者Anti Earth