中等网格分辨率下Gmsh细长体网格划分失败问题

问题说明 #

在开展Gmsh网格收敛性研究时，当网格细化到某一级别（网格尺寸仍远大于所有可配置公差，本不应触发浮点/几何公差问题），2D表面网格划分阶段会弹出如下错误：

Error : Identical points in triangulation: increase element size or Mesh.RandomFactor

尝试过增大单元尺寸（但收敛性研究需要低分辨率，此方法不可行）、修改Mesh.RandomFactor、调整Geometry.Tolerance、Mesh.ToleranceInitialDelaunay或Geometry.ReparamOnFaceRobust参数，均无法解决问题。

可复现代码如下（当前SMIN和SMAX的设置会触发错误，增大尺寸则可消除）：

import numpy as np
import gmsh

R_INNER = 0.001
R_OUTER = 0.003
LENGTH  = 0.500
ANGLE   = 2.0          # 楔形角度，单位：度
SMIN    = 1e-5
SMAX    = 5e-5

gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 5)
gmsh.model.add("cylinder_wedge")
occ = gmsh.model.occ

# x-z平面内的横截面（x为轴向，z为径向，y=0）
# 绕x轴旋转生成薄楔形
p = [
    occ.add_point(0,      0, R_INNER),   
    occ.add_point(LENGTH, 0, R_INNER),  
    occ.add_point(LENGTH, 0, R_OUTER),   
    occ.add_point(0,      0, R_OUTER),   
]
lines = [
    occ.add_line(p[0], p[1]),   
    occ.add_line(p[1], p[2]),   
    occ.add_line(p[2], p[3]), 
    occ.add_line(p[3], p[0]),
]
face = occ.add_plane_surface([occ.add_curve_loop(lines)])

revolved = occ.revolve(
    [(2, face)],
    0, 0, 0,            
    1, 0, 0,            
    np.deg2rad(ANGLE),
)
occ.synchronize()

vol = next(tag for dim, tag in revolved if dim == 3)
gmsh.model.add_physical_group(3, [vol], name="VOLUME")

gmsh.option.setNumber("Mesh.MeshSizeMin",             SMIN)
gmsh.option.setNumber("Mesh.MeshSizeMax",             SMAX)
gmsh.option.setNumber("Mesh.Algorithm",               6)   # Frontal-Delaunay
gmsh.option.setNumber("Geometry.ReparamOnFaceRobust", 1)

gmsh.model.mesh.generate(3)

gmsh.write("cylinder.msh")
gmsh.finalize()

原因分析 #

这个错误的核心是：旋转生成的薄楔形侧面（仅2度夹角）在网格细化时，由于几何特征过于扁平，Frontal-Delaunay三角化算法在离散过程中出现了数值上近乎重合的点，触发了Gmsh的几何检查机制。

可行解决方案 #

1. 为薄侧面指定专用网格算法 #

强制为旋转生成的小角度侧面使用更适合狭长/扁平几何的网格算法（比如MeshAdapt），替代默认的Frontal-Delaunay：
在occ.synchronize()之后添加以下代码：

# 获取所有2D面
faces = [tag for dim, tag in gmsh.model.getEntities(2)]
for face_tag in faces:
    # 为每个面设置MeshAdapt算法（算法编号1）
    gmsh.model.mesh.setAlgorithm(2, face_tag, 1)

2. 局部控制网格尺寸 #

不为整个模型设置统一的最小/最大尺寸，而是针对不同区域单独配置，避免薄侧面被过度细化：
替换全局尺寸设置的代码，改为局部尺寸控制：

# 为径向的两条线设置最小网格尺寸
radial_lines = [lines[1], lines[3]]
for line_tag in radial_lines:
    gmsh.model.mesh.setSize([(1, line_tag)], SMIN)

# 为轴向的两条线设置最大网格尺寸
axial_lines = [lines[0], lines[2]]
for line_tag in axial_lines:
    gmsh.model.mesh.setSize([(1, line_tag)], SMAX)

# 保留其他网格选项
gmsh.option.setNumber("Mesh.Algorithm",               6)   # Frontal-Delaunay
gmsh.option.setNumber("Geometry.ReparamOnFaceRobust", 1)

3. 清理几何重合特征 #

在同步几何模型前，显式清理可能存在的重复几何实体：
在occ.synchronize()之前添加：

occ.removeAllDuplicates()

修改后的完整代码示例 #

整合上述方案后的代码如下：

import numpy as np
import gmsh

R_INNER = 0.001
R_OUTER = 0.003
LENGTH  = 0.500
ANGLE   = 2.0          # 楔形角度，单位：度
SMIN    = 1e-5
SMAX    = 5e-5

gmsh.initialize()
gmsh.option.setNumber("General.Verbosity", 5)
gmsh.model.add("cylinder_wedge")
occ = gmsh.model.occ

# x-z平面内的横截面（x为轴向，z为径向，y=0）
# 绕x轴旋转生成薄楔形
p = [
    occ.add_point(0,      0, R_INNER),   
    occ.add_point(LENGTH, 0, R_INNER),  
    occ.add_point(LENGTH, 0, R_OUTER),   
    occ.add_point(0,      0, R_OUTER),   
]
lines = [
    occ.add_line(p[0], p[1]),   
    occ.add_line(p[1], p[2]),   
    occ.add_line(p[2], p[3]), 
    occ.add_line(p[3], p[0]),
]
face = occ.add_plane_surface([occ.add_curve_loop(lines)])

revolved = occ.revolve(
    [(2, face)],
    0, 0, 0,            
    1, 0, 0,            
    np.deg2rad(ANGLE),
)

# 清理重复几何实体
occ.removeAllDuplicates()
occ.synchronize()

vol = next(tag for dim, tag in revolved if dim == 3)
gmsh.model.add_physical_group(3, [vol], name="VOLUME")

# 局部控制网格尺寸
# 为径向的两条线设置最小网格尺寸
radial_lines = [lines[1], lines[3]]
for line_tag in radial_lines:
    gmsh.model.mesh.setSize([(1, line_tag)], SMIN)

# 为轴向的两条线设置最大网格尺寸
axial_lines = [lines[0], lines[2]]
for line_tag in axial_lines:
    gmsh.model.mesh.setSize([(1, line_tag)], SMAX)

# 为薄侧面指定MeshAdapt算法
faces = [tag for dim, tag in gmsh.model.getEntities(2)]
for face_tag in faces:
    gmsh.model.mesh.setAlgorithm(2, face_tag, 1)

gmsh.option.setNumber("Geometry.ReparamOnFaceRobust", 1)

gmsh.model.mesh.generate(3)

gmsh.write("cylinder.msh")
gmsh.finalize()

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