如何用Python基于顶点坐标生成棱柱侧面四边形非结构化表面网格

我有一个如下图所示的长方体棱柱，已知其所有顶点的坐标。

我希望将该棱柱的侧面离散为采用指定目标边长四边形单元的非结构化表面网格。具体需求为：无需将网格保存为文件，仅需得到表示网格单元的列表，列表中每个单元需包含构成该单元的节点引用（同时涵盖节点ID与节点坐标信息），请问是否有Python库可实现该功能？

对于上图所示的规则长方体棱柱，仅通过调整坐标结合numpy.linspace即可非常简便地生成结构化表面网格。但我还需要支持对带锥度、或上下底面存在差异的棱柱进行网格划分，这类场景下结构化网格容易产生畸变程度很高的单元，而我希望生成的单元尽可能接近边长等于目标长度的正方形。

理想情况下，Python脚本首先定义待划分网格的顶点坐标以及需要划分的面，代码如下：

import numpy as np
# Geometric properties
aspect_ratio = 9.
c = 1e3 # [mm]
b = aspect_ratio*c
non_dimensional_height = 0.1
h = c*non_dimensional_height # [mm]
non_dimensional_thickness = 0.001
t = c*non_dimensional_thickness # [mm]
# Objective length of the side of the quadrilateral elements
side_length = 10 # [mm]
# Define points coordinates
points_xyz = np.array([(0., 0., h/2),   # 0 -> A
                       (c, 0., h/2),    # 1 -> B
                       (0., 0., -h/2),  # 2 -> C
                       (c, 0., -h/2),   # 3 -> D
                       (0., b, h/2),    # 4 -> A'
                       (c, b, h/2),     # 5 -> B'
                       (0., b, -h/2),   # 6 -> C'
                       (c, b, -h/2)])   # 7 -> D'
# Define faces by sequence of points
faces = np.array([(0,1,5,4), (1,2,6,5), (2,3,7,6), (3,0,4,7)])

后续逻辑需完成各个面的离散，返回包含所有非结构化网格单元的列表，每个单元关联其对应的节点信息（同时包含节点ID与坐标）。

我了解到pyvista、gmsh等相关库，但印象中这类库更偏向于处理从文件导入的网格，而非直接基于输入坐标生成网格、输出网格连通性信息。

回答 #

你对gmsh和pyvista的认知存在偏差，这两个库完全支持直接从内存中的坐标数据生成网格、直接读取连通性和节点信息，不需要经过文件读写。针对四边形主导非结构化表面网格的需求，最简便的方案是用gmsh Python API直接构造平面面域，设置目标单元尺寸后生成网格，再提取节点和单元数据即可。

具体实现逻辑：

1. 初始化gmsh实例，关闭冗余日志输出
2. 逐点添加所有顶点坐标，记录每个顶点对应的gmsh点标签
3. 逐面根据顶点标签构造四边形面的闭合线环，生成平面曲面
4. 统一设置网格单元尺寸为指定目标边长，选用四边形优先的表面网格算法
5. 生成2D表面网格后，直接从内存读取节点坐标、单元连通性数据，整理成需要的格式
6. 清空gmsh实例，全程不需要写入任何文件

对应给出的起始代码，可直接运行的完整实现如下：

import numpy as np
import gmsh

# ---------------------- 原有参数定义 ----------------------
aspect_ratio = 9.
c = 1e3 # [mm]
b = aspect_ratio*c
non_dimensional_height = 0.1
h = c*non_dimensional_height # [mm]
non_dimensional_thickness = 0.001
t = c*non_dimensional_thickness # [mm]
side_length = 10 # [mm]
points_xyz = np.array([(0., 0., h/2),   # 0 -> A
                       (c, 0., h/2),    # 1 -> B
                       (0., 0., -h/2),  # 2 -> C
                       (c, 0., -h/2),   # 3 -> D
                       (0., b, h/2),    # 4 -> A'
                       (c, b, h/2),     # 5 -> B'
                       (0., b, -h/2),   # 6 -> C'
                       (c, b, -h/2)])   # 7 -> D'
# 注意：使用时需保证每个面的顶点按顺时针/逆时针连续闭合，避免生成扭曲面
faces = np.array([(0,1,3,2), (1,5,7,3), (5,4,6,7), (4,0,2,6)]) # 修正为长方体四个侧面的正确点序
# -----------------------------------------------------------

# 初始化gmsh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # 关闭控制台冗余日志
model = gmsh.model
model.add("prism_mesh")

# 添加所有顶点
point_tags = []
for x,y,z in points_xyz:
    tag = model.geo.addPoint(x, y, z, meshSize=side_length)
    point_tags.append(tag)

# 添加所有待划分面
for face_pid in faces:
    # 逐边构造线
    line_tags = []
    for i in range(4):
        p1 = point_tags[face_pid[i]]
        p2 = point_tags[face_pid[(i+1)%4]]
        line_tags.append(model.geo.addLine(p1, p2))
    # 构造闭合线环与平面曲面
    loop_tag = model.geo.addCurveLoop(line_tags)
    model.geo.addPlaneSurface([loop_tag])

model.geo.synchronize()

# 设置网格生成参数：生成四边形主导网格，控制单元尺寸尽可能接近目标值
gmsh.option.setNumber("Mesh.Algorithm", 8) # 选用Delaunay四边形网格算法，单元正方形占比最高
gmsh.option.setNumber("Mesh.RecombineAll", 1) # 所有面自动重组为四边形单元
gmsh.option.setNumber("Mesh.MeshSizeMin", side_length*0.9)
gmsh.option.setNumber("Mesh.MeshSizeMax", side_length*1.1) # 单元尺寸波动控制在10%以内

# 生成2D表面网格
model.mesh.generate(2)

# 提取网格数据整理为目标格式
# 提取所有节点，生成ID到坐标的映射
node_tags, node_coords, _ = model.mesh.getNodes()
node_coords = node_coords.reshape(-1, 3)
node_id_map = {tag:i for i, tag in enumerate(node_tags)}
node_list = [{"node_id": i, "coord": coord} for i, coord in enumerate(node_coords)]

# 提取所有四边形单元
element_list = []
elem_id = 0
for surf_dim, _ in model.getEntities(2):
    elem_types, _, elem_node_tags = model.mesh.getElements(surf_dim)
    for etype, enodes in zip(elem_types, elem_node_tags):
        if etype == 3: # etype=3对应4节点四边形单元
            enodes_arr = np.array(enodes).reshape(-1,4)
            for cell_nodes in enodes_arr:
                cell_node_ids = [node_id_map[n] for n in cell_nodes]
                element_list.append({
                    "element_id": elem_id,
                    "nodes": [node_list[nid] for nid in cell_node_ids]
                })
                elem_id +=1

# 清空gmsh实例，无任何文件读写操作
gmsh.finalize()

# 结果验证
print(f"生成节点总数：{len(node_list)}")
print(f"生成四边形单元总数：{len(element_list)}")

使用说明：

• 自定义棱柱顶点与面时，仅需保证每个面的顶点共面、按顺时针/逆时针顺序连续闭合即可，锥度、异形面都可以正常生成网格
• 可通过调整Mesh.MeshSizeMin和Mesh.MeshSizeMax的系数放宽/收紧单元尺寸的波动范围
• 如果需要更高的四边形质量，可设置Mesh.RecombineMinimumQuality参数（取值0~1），低于质量阈值的单元会自动退化为三角形避免畸变
• 最终输出的element_list就是需要的单元列表，每个单元同时包含节点ID与三维坐标信息

如果不想依赖gmsh的几何内核，也可以用triangle库配合四边形重组逻辑实现，但gmsh对异形棱柱的网格平滑、质量优化更成熟，不需要额外开发单元校验相关逻辑。

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