Gmsh Python API中getElementsByType无法返回全部单元求助

在Python中使用Gmsh的getElementsByType函数获取矩形域（t1.py示例）的所有三角形单元时，网格实际有796个三角形，但该函数仅返回726个，标签1到70的单元未被获取到。相关代码如下：

import gmsh
import sys
import numpy as np

gmsh.initialize()

gmsh.model.add("t1")
lc = 1e-2

gmsh.model.geo.addPoint(0, 0, 0, lc, 1)
gmsh.model.geo.addPoint(.1, 0, 0, lc, 2)
gmsh.model.geo.addPoint(.1, .3, 0, lc, 3)
p4 = gmsh.model.geo.addPoint(0, .3, 0, lc)

gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(3, 2, 2)
gmsh.model.geo.addLine(3, p4, 3)
gmsh.model.geo.addLine(4, 1, p4)

gmsh.model.geo.addCurveLoop([4, 1, -2, 3], 1)

gmsh.model.geo.addPlaneSurface([1], 1)

gmsh.model.geo.synchronize()

gmsh.model.addPhysicalGroup(1, [1, 2, 4], 5)
gmsh.model.addPhysicalGroup(2, [1], name="My surface")

gmsh.model.mesh.generate(2)

gmsh.write("t1.msh")
gmsh.open("t1.msh")

dim = -1
tag = -1
# we print the coordinates of nodes
nodeTags, coords, paraCoord = gmsh.model.mesh.getNodes(dim, tag)
coords = coords.reshape(-1,3)
xy = coords.reshape(-1,3)
for Tag, xy_e in zip(nodeTags, xy):
    print(f"Node #{Tag} is at {xy_e}")

defined_element_type = gmsh.model.mesh.getElementTypes()
print(defined_element_type)

eletype = 2
tag = -1
eleTags, nodeTags = gmsh.model.mesh.getElementsByType(eletype, tag)
nodeTags = nodeTags.reshape(-1,3)

nodesEachEle = []
for Tag, nodes in zip(eleTags, nodeTags):
    print(f"Element #{Tag} has nodes {nodes}")
    nodesEachEle.append(nodes)
#if '-nopopup' not in sys.argv:
    #gmsh.fltk.run()

gmsh.finalize()

问题原因 #

核心问题出在生成网格后调用的gmsh.write和gmsh.open操作：

• 当执行gmsh.open("t1.msh")时，Gmsh会创建新的模型实例加载网格文件，直接替换掉之前生成网格的原始模型上下文。
• 新加载的模型中，单元标签体系与原始生成的模型不一致，标签1-70的三角形单元属于边界关联的网格单元，在新模型中未被纳入getElementsByType的查询范围，导致丢失。

解决方案 #

直接删除gmsh.write("t1.msh")和gmsh.open("t1.msh")两行代码，生成网格后直接在当前模型上下文内查询单元即可，无需写入后重新打开文件。

修改后的核心代码片段：

gmsh.model.mesh.generate(2)

# 移除write和open操作，直接在当前模型中查询单元
dim = -1
tag = -1
nodeTags, coords, paraCoord = gmsh.model.mesh.getNodes(dim, tag)
coords = coords.reshape(-1,3)
xy = coords.reshape(-1,3)
for Tag, xy_e in zip(nodeTags, xy):
    print(f"Node #{Tag} is at {xy_e}")

# 后续查询代码保持不变

补充说明：如果需要保存网格文件，仅保留gmsh.write("t1.msh")即可，无需重新打开——重新打开会重置模型状态，破坏原始生成的单元标签与结构关联。

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