Python开源网格工具 | Delaunay网格生成 | 开源库 + Bowyer-Watson算法

Delaunay 三角剖分是计算几何中的经典算法，广泛应用于网格生成、地形建模和计算机图形学。本文将深入解析 PyMeshGen 中的 Delaunay 模块，揭示如何基于 Gmsh 的 Bowyer-Watson 算法实现高质量的二维三角网格生成。

一、引言

在计算流体力学（CFD）和有限元分析（FEA）中，三角网格是最基础的非结构网格类型。Delaunay 三角剖分因其独特的空圆特性（Empty Circumcircle Property）而成为网格生成的首选方法：

• 最优角度：最大化最小内角，避免狭长三角形
• 唯一性：给定一组点，Delaunay 三角剖分在非退化情况下唯一
• 局部优化：任意两三角形组成的四边形中，对角线是最优选择

PyMeshGen 的 delaunay\ 模块实现了完整的二维 Delaunay 三角剖分系统，支持两种后端：

1. Bowyer-Watson 后端：基于 Gmsh 风格的增量插入算法
2. Triangle 后端：封装 Jonathan Shewchuk 的经典 Triangle 程序

本文将重点解析 Bowyer-Watson 后端的实现原理。

注：网格质量仍有优化空间，但是边界恢复正确。

二、算法总体框架

2.1 模块架构

PyMeshGen 的 Delaunay 模块采用分层设计，将输入归一化、核心算法、后处理解耦：

delaunay/├── __init__.py              # 包入口，公开接口├── bw_utils.py              # 输入归一化与后端分发├── bw_core_stable.py        # Bowyer-Watson 主算法├── bw_cavity.py             # 空腔搜索与重连├── bw_types.py              # Gmsh 风格数据结构├── bw_predicates.py         # 几何谓词├── triangle_backend.py      # Triangle 后端├── postprocess.py           # 轻量后处理└── validation.py            # 测试验证工具

分层职责：

2.2 典型调用链

Delaunay 模块并不是独立入口，而是 core.py 的一个子系统：

core.generate_mesh()  -> QuadtreeSizing(...)                    # 构建尺寸场  -> create_bowyer_watson_mesh(...)         # delaunay 公共入口      -> bw_utils._build_boundary_input()   # front -> 点/边/孔洞      -> backend dispatch         -> BowyerWatsonMeshGenerator       # Bowyer-Watson 后端         -> create_triangle_mesh            # Triangle 后端  -> core._recover_delaunay_boundary_edges()  # 轻量边翻转补恢复  -> Unstructured_Grid.from_cells()           # 统一转网格对象

2.3 公共入口

业务层只需调用一个函数：

from delaunay import create_bowyer_watson_meshpoints, simplices, boundary_mask = create_bowyer_watson_mesh(    boundary_front=front,    sizing_system=sizing,    backend="bowyer_watson",  # 或 "triangle")

该入口负责：

1. 接收 boundary_front 和 QuadtreeSizing
2. 归一化边界输入为 BoundaryInput
3. 自动识别孔洞与外边界
4. 根据 backend 参数选择后端
5. 统一返回 (points, simplices, boundary_mask) 三元组

设计原则：上层不需要知道 Bowyer-Watson 内部拓扑如何组织，只关心统一数组输出；后端切换不改变返回格式。

三、边界输入归一化

3.1 BoundaryInput 结构

算法核心不直接理解 front 对象，而是通过 BoundaryInput 接收归一化的边界信息：

BoundaryInput├── boundary_points   # 去重后的边界点坐标├── boundary_edges    # 边界边索引对├── holes             # 孔洞环列表└── outer_boundary    # 主外边界点环

3.2 归一化流程

boundary_front  -> _build_front_graph()           # 以 node.hash 为图节点，阵面为无向边  -> _trace_boundary_loops()        # 从 front 图中提取闭合 loop  -> _classify_boundary_loops()     # 基于面积与重心包含关系判别外边界/孔洞  -> _extract_boundary_points_and_edges()  # 索引稳定化  -> BoundaryInput

关键处理点：

1. front 图重建：以 node.hash 为图节点，以阵面为无向边构建图结构
2. 边界环追踪：从 front 图中提取所有闭合的边界环
3. 外边界/孔洞判别：基于面积大小与重心包含关系，面积最大的环为外边界，被外边界包含的环为孔洞
4. 索引稳定化：所有边界边在进入核心算法前重映射到连续索引

这样做的好处是：算法核心不需要理解 front 对象；Triangle 后端与 Bowyer-Watson 后端可共享同一边界输入；测试与调试时可以直接输出边界点/边，不依赖原始 front 结构。

四、核心数据结构

4.1 MTri3 —— Gmsh 风格的三角形包装类

MTri3 是 Bowyer-Watson 主路径的核心数据结构，直接参考 Gmsh 的 MTri3 类设计：

classMTri3:    __slots__ = ['vertices',        # tuple: 排序后的顶点索引 (v0, v1, v2)'neighbors',       # list[MTri3|None]: 三个邻居三角形'circumcenter',    # np.ndarray: 外接圆心 [x, y]'circumradius',    # float: 外接圆半径'deleted',         # bool: 懒删除标记'idx',             # int: 三角形唯一标识'quality',         # float: 质量度量 (0-1)    ]

关键设计：

1. 懒删除（Lazy Deletion）

• 空腔搜索、边翻转、局部修复时避免频繁从容器中物理移除
• 通过 deleted 标记逻辑删除，统一压缩时再物理移除

2. 显式邻接（Explicit Adjacency）

• neighbors[3] 维护三角形邻接关系，支持 O(1) 的邻居访问
• Cavity 搜索和边恢复可以直接走邻接关系，无需重建边索引

3. 缓存几何量

• circumcenter 和 circumradius 缓存外接圆信息
• 避免重复计算，circumradius 同时作为三角形质量度量

4. 顶点索引排序

• vertices 自动按升序存储，简化三角形相等性判断
• __eq__ 和 __hash__ 基于顶点集合

4.2 EdgeXFace —— 边-面关系结构

EdgeXFace 表示"边-面"关系，主要用于空腔 shell 收集和约束边标识：

classEdgeXFace:    __slots__ = ['triangle', 'local_edge_idx', 'is_constrained', 'neighbor']

• triangle + local_edge_idx 唯一标识一条边
• is_constrained 标记约束边，防止被删除或翻转
• neighbor 指向相邻的 EdgeXFace

4.3 TriangulationState —— 常驻拓扑索引层

TriangulationState 是在 MTri3 之上的常驻索引缓存，用于替代反复现建现用的临时 edge map / vertex map：

classTriangulationState:    __slots__ = ['triangles',           # 三角形列表'triangle_by_id',      # ID -> 三角形映射'edge_to_tris',        # 边 -> 关联三角形列表'vertex_to_tris',      # 顶点 -> 关联三角形列表'vertex_to_neighbors', # 顶点 -> 邻点集合    ]

职责：

1. 邻接重建：用一次边遍历同时完成 neighbors 回填和 edge index 构建
2. 高频查询：查询边是否存在、某边 incident triangles、某顶点 incident triangles、某顶点邻点集合
3. 懒删除压缩：在需要时统一物理压缩 deleted triangles，并刷新索引

几乎所有局部操作都依赖这层索引：cavity 扩张、边翻转、边界恢复、局部 strip 搜索。

五、几何谓词与数值稳健性

5.1 核心谓词

Delaunay 三角剖分依赖两个基本几何谓词：

1. orient2d —— 方向测试

判断点 p 相对于有向直线 (a→b) 的位置：

deforient2d(ax, ay, bx, by, px, py) -> float:"""返回：      > 0：p 在直线左侧      = 0：p 在直线上      < 0：p 在直线右侧    """return (bx - ax) * (py - ay) - (by - ay) * (px - ax)

2. incircle —— 圆内测试

判断点 p 是否在三角形 (a, b, c) 的外接圆内：

defincircle(ax, ay, bx, by, cx, cy, px, py) -> float:"""返回：      > 0：p 在圆内      = 0：p 在圆上      < 0：p 在圆外    """    adx = ax - px; ady = ay - py    bdx = bx - px; bdy = by - py    cdx = cx - px; cdy = cy - py    alift = adx*adx + ady*ady    blift = bdx*bdx + bdy*bdy    clift = cdx*cdx + cdy*cdy    det = (adx*(bdy*clift - cdy*blift)          - ady*(bdx*clift - cdx*blift)          + alift*(bdx*cdy - cdx*bdy))    orient = (bx-ax)*(cy-ay) - (by-ay)*(cx-ax)return det * orient

5.2 稳健性策略

几何谓词的数值稳定性是 Delaunay 算法的核心挑战。PyMeshGen 采用多层防护：

为什么单独拆文件：

1. 算法逻辑与几何数值分离
2. 可以在不同实现中复用
3. 便于将来替换更强的鲁棒谓词实现

六、Bowyer-Watson 核心算法

6.1 算法概述

Bowyer-Watson 算法是一种增量插入的 Delaunay 三角剖分方法，核心思想是：

1. 创建超级三角形（包含所有点）2. 逐个插入新点：   a. 找到包含新点的三角形   b. 搜索空腔（所有外接圆包含新点的三角形）   c. 删除空腔内的三角形   d. 用空腔边界与新点重新连接3. 删除含超级顶点的三角形

6.2 BowyerWatsonMeshGenerator 主流程

BowyerWatsonMeshGenerator 是 backend="bowyer_watson" 的唯一实现，设计上更接近 Gmsh meshGFaceDelaunayInsertion.cpp 的二维思路。

主入口阶段顺序：

阶段 1: 初始三角剖分  ├─ 创建超级三角形  ├─ 按顺序插入所有边界点  ├─ 使用 cavity 搜索删除被新点破坏的三角形  ├─ 用 shell 边与新点重新连接  ├─ 删除含超级顶点的三角形  └─ 立即恢复初始剖分中缺失的受保护边阶段 2: Gmsh 风格迭代插点  ├─ 优先级队列（按外接圆半径排序）  ├─ 尺寸场驱动（target_size = sizing_system.spacing_at(tri_center)）  ├─ 质量阈值（边界区与内区使用不同阈值）  ├─ 早停策略（定期全量检查剩余不满足要求的单元数量）  └─ 压缩与重建（deleted 积累到一定规模后压缩容器并重建队列）阶段 2.5: CDT 边界恢复  └─ swap-first 的精确受保护边恢复阶段 2.6~2.95: 孔洞/域外清理 + 再恢复  ├─ 清理孔洞内三角形  ├─ 清理域外三角形  └─ 必要时再次做 CDT 恢复阶段 3: Laplacian 平滑  └─ 节点位置优化阶段 3.5: 重叠/重复/退化/交叉清理  ├─ 重叠三角形清理  ├─ 重复三角形清理  ├─ 严格相交三角形清理  ├─ 退化三角形清理  ├─ 被隔离边界点修复  └─ 孔洞侧 boundary-fan 清理输出压缩与导出

6.3 初始三角剖分

初始剖分的目标不是得到最终高质量网格，而是建立一个：

• 覆盖整个域
• 邻接正确
• 边界点完整
• 可继续细化的初始三角网格

# 伪代码def_initial_triangulation(self):# 1. 创建超级三角形    super_tri = self._create_super_triangle(boundary_points)# 2. 逐个插入边界点for point in boundary_points:        cavity = self._find_cavity(point)        shell = self._collect_cavity_shell(cavity)        self._delete_cavity(cavity)        self._reconnect(shell, point)# 3. 删除含超级顶点的三角形    self._remove_super_triangles()# 4. 恢复缺失的受保护边    self._recover_protected_edges()

6.4 迭代细化策略

细化主循环 _insert_points_iteratively() 的核心元素：

1. 优先级队列

# 按外接圆半径排序，优先处理"坏三角形"bad_triangles = sorted(triangles, key=lambda t: t.circumradius, reverse=True)

2. 尺寸场驱动

target_size = sizing_system.spacing_at(tri_center)if tri.circumradius > target_size:# 需要细化

3. 质量阈值

边界区与内区使用不同阈值：

• 边界区：更宽松，减少边界恢复压力
• 内区：更严格，保证网格质量

4. 早停策略

定期全量检查剩余不满足要求的单元数量，避免无限循环。

5. 压缩与重建

deleted 三角形积累到一定规模后，压缩容器并重建队列，避免内存膨胀。

6.5 细化点选择

候选点选择由 _select_refinement_point() 完成：

def_select_refinement_point(self, triangle):# 1. 若是靠近受保护边的低质量三角形，优先尝试 off-centerif self._near_protected_edge(triangle):        candidate = self._off_center(triangle)# 2. 否则默认使用 circumcenterelse:        candidate = triangle.circumcenter# 3. 若候选点不满足域内/孔洞/最小间距要求，则拒绝ifnot self._is_valid_insertion(candidate):returnNonereturn candidate

off-center 策略：对边界附近使用更保守的插点，减少边界恢复压力；对内部维持 Delaunay 细化效率。

七、空腔搜索与重连

7.1 Cavity 搜索算法

bw_cavity.py 实现空腔算法：

start_tri  -> recur_find_cavity()      -> 沿邻接递归扩张      -> 碰到受保护边时停止      -> 收集 cavity triangles      -> 收集 shell edges  -> insert_vertex()      -> 新点连接 shell      -> 建立新三角形      -> 更新邻接

核心约束：

• 受保护边不能跨越
• shell 必须能包围 cavity
• 失败时必须恢复 deleted 标记和临时点

7.2 Shell 收集

空腔边界（Shell）的收集是重连的关键：

defcollect_cavity_shell(cavity_triangles: List[MTri3]) -> List[EdgeXFace]:"""收集空腔的边界边（Shell）。    算法：    - 遍历空腔中所有三角形的所有边    - 如果某条边只被一个空腔三角形拥有，则是边界边    """    edge_count = {}    edge_to_tris = {}for tri in cavity_triangles:for i in range(3):            edge_key = tri.get_edge_sorted(i)            edge_count[edge_key] = edge_count.get(edge_key, 0) + 1            edge_to_tris[edge_key] = (tri, i)# 边界边是只出现一次的边    shell = []for edge_key, count in edge_count.items():if count == 1:            tri, local_idx = edge_to_tris[edge_key]            shell.append(EdgeXFace(tri, local_idx))return shell

7.3 失败回滚

_rollback_failed_cavity_insertion() 与 _insert_refinement_point() 负责：

• 回滚 deleted 三角形
• 删除失败插入产生的临时点
• 防止留下"孤儿点""空洞"或不完整局部拓扑

这是当前实现区别于早期版本的重要稳健性增强点。

八、边界恢复

8.1 两层边界恢复

模块中存在两层边界恢复：

1. 核心恢复：bw_core_stable.py

def_constrained_delaunay_triangulation(self):"""面向受保护边执行 swap-first 的精确恢复"""for protected_edge in self.protected_edges:ifnot self._edge_exists(protected_edge):            self._recover_edge(protected_edge)

2. 轻量恢复：postprocess.py

defrecover_boundary_edges_by_swaps(points, simplices, boundary_edges):"""对数组级三角形做边翻转恢复"""# 在 core.py 中作为数组级补恢复使用

设计目的：核心算法内部尽量保证精确约束边；上层在最终三角数组层面还有一次便宜的补救机会。

8.2 边界恢复顺序

Gmsh 主路径中的孔洞处理顺序非常关键：

1. 先完整细化2. 再恢复边界3. 再清理孔洞与域外三角形4. 必要时再次做 CDT 恢复

原因：

• 若过早删洞，局部 cavity 不完整，边界恢复会更困难
• 若清理后不再恢复，某些真实边界边可能被再次丢失

九、孔洞与域清理

9.1 孔洞识别与清理

孔洞清理的核心是判断三角形是否在孔洞内部：

def_is_in_hole(self, triangle, holes):"""判断三角形重心是否在任意孔洞内"""    center = triangle.circumcenterfor hole in holes:if self._point_in_polygon(center, hole):returnTruereturnFalse

9.2 多边界环场景

当前实现已经显式支持通过边界连通分量区分不同边界环：

• 对同一边界环上的 boundary fan 三角形可放宽尺寸细化
• 对跨外边界/孔洞边界的桥接三角形仍允许内部细化

这一点对环域类算例（如 quad_quad）尤其关键。

十、拓扑清理

10.1 最终输出前的清理

最终输出前，Gmsh 主路径会执行：

def_cleanup_topology(self):# 1. 重叠三角形清理    self._remove_overlapping_triangles()# 2. 重复三角形清理    self._remove_duplicate_triangles()# 3. 严格相交三角形清理    self._remove_intersecting_triangles()# 4. 退化三角形清理    self._remove_degenerate_triangles()# 5. 被隔离边界点修复    self._repair_isolated_boundary_points()# 6. 孔洞侧 boundary-fan 清理    self._cleanup_boundary_fans()

这一步是"结果正确性"的最后保障层。

10.2 拓扑验证

validation.py 提供单元测试和验收检查：

defcheck_boundary_edges(points, simplices, boundary_edges):"""检查边界边是否完全恢复"""defcheck_hole_cleanup(points, simplices, holes):"""检查孔洞是否完全清理"""defcheck_topology_clean(points, simplices):"""检查是否存在严格交叉"""

核心思路：重新解析输入 CAS 边界，将输入边界与输出网格建立节点映射，检查边界是否恢复、孔洞是否为空、是否存在严格交叉。

十一、Triangle 后端

11.1 设计目标

triangle_backend.py 为 Jonathan Shewchuk 的 Triangle 提供本地包装：

• 不修改第三方源码
• 使用 CLI 而非 DLL/ctypes 桥接
• 继续复用项目已有尺寸场和边界 front 体系

11.2 工作流

QuadtreeSizing  -> 内部点采样  -> 写 mesh.poly  -> 调用 triangle.exe  -> 解析 mesh.1.node / mesh.1.ele  -> 返回 (points, simplices, boundary_mask)

11.3 内部点采样策略

支持两种策略：

Triangle 后端本质上是：边界与尺寸场仍由 PyMeshGen 控制，三角剖分核心委托给 Triangle。

十二、后处理与验证

12.1 postprocess.py

提供两类轻量工具：

1. recover_boundary_edges_by_swaps()

对数组级三角形做边翻转恢复，不依赖 MTri3。

2. is_topology_valid()

检查：

• 是否存在超过 2 个单元共边
• 网格连通性
• 严格边交叉

它的定位是：便于 core.py 在统一结果层面做补检查。

12.2 与主流程的集成

core.py 里 mesh_type=4 的职责分工如下：

# core.py 伪代码defgenerate_mesh():# 1. 读取输入网格并构造初始 front    front = self._build_initial_front()# 2. 构造 QuadtreeSizing    sizing = QuadtreeSizing(front, params)# 3. 决定是否先生成边界层if has_boundary_layer:        front = self._generate_boundary_layer(front)# 4. 调用 create_bowyer_watson_mesh()    points, simplices, boundary_mask = create_bowyer_watson_mesh(        boundary_front=front,        sizing_system=sizing,        backend=params.delaunay_backend,    )# 5. 对非 Triangle 后端做轻量边翻转恢复if params.delaunay_backend != "triangle":        points, simplices = recover_boundary_edges_by_swaps(...)# 6. 转成 Unstructured_Grid    mesh = Unstructured_Grid.from_cells(points, simplices)return mesh

十三、配置参数与行为控制

与 delaunay\ 强相关的参数主要有：

行为补充规则：带边界层时，上层可能强制切换到 Triangle 后端。

十四、总结

PyMeshGen 的 Delaunay 模块的核心设计思想可以概括为：

1. 输入统一：front 先归一化成点/边/孔洞
2. 后端可切换：Bowyer-Watson 与 Triangle 共用同一上层接口
3. 核心算法分层：几何谓词、cavity、数据结构、后处理分开实现
4. 结果优先：边界恢复、孔洞清理、拓扑清理构成多层保障
5. 集成友好：始终以 (points, simplices, boundary_mask) 作为对上层的稳定契约

这套设计使得模块既能服务于当前项目的工程化网格生成，又保留了继续对齐 Gmsh / Triangle / 本项目自定义策略的扩展空间。