浅谈曲面网格生成的理论基础

taotao_2016 2022-03-28

展开全文

俄乌战争的爆发深刻地改变了历史的进程。可以想见，未来计算机科学将会更加侧重服务于实体工业。在计算机辅助设计（CAD），计算机辅助工程（CAE）和计算机辅助制造（CAM）领域中，工业软件起到了决定性的作用。工业软件的底层逻辑如下：人们在各种工业产品的设计过程中需要对其模拟仿真，以便发现缺陷修改设计方案。模拟仿真核心是用计算机来求解各种物理微分、积分方程，来预测机械的物理行为，例如电磁场的分布变化、热力学现象、应力应变、流体现象等等。在具有复杂拓扑和几何的曲面或者实体上求解偏微分方程，通常采用有限元法、有限体积法或者等几何分析方法，这些方法都强烈依赖于曲面实体的网格生成。有限元法、有限体积法需要非结构化的三角剖分，也可以基于结构化网格生成（例如四边形、六面体剖分）以进一步提高计算精度。而等几何分析方法完全依赖于结构化网格生成。由此，通常的工业软件有三个主要功能模块：网格生成、微分方程求解和可视化。在实际应用中，网格生成占据了73%左右的时间成本和经济成本。由此可见，网格生成是现代工业软件的核心基础。

虽然网格生成技术至关重要，但是相比于数值偏微分方程求解和可视化技术，其成熟程度远远低于其他两者，很多非常基本的理论问题依然长期悬而未决。在所有计算机科学、应用数学、机械力学等院系中，数学物理方程、数值代数、数值微分方程都是基础课程，计算机图形学和可视化也是计算机系的必修课。经典的数值方法都已经被Matlab等软件所容纳，变成了一条条简单指令；传统的图形渲染算法也固化到GPU的硬件里面，变成简单函数调用。但是，很少有大学开设网格生成课程，即便有类似研究生课程，也没有教材可以囊括完整的知识体系，而是东鳞西爪，零星算法的集成。虽然工业界有巨大需求，但是因为这一领域具有悠久的历史，容易解决的问题早已解决，余下的问题过于艰深，真正的突破需要比较深刻的现代拓扑几何理论。无论是学习理论知识还是工程积累，都需要板凳一坐十年冷。对于计算机科学领域的年轻学者而言，贸然投入研究，风险过大，因此这一领域在学术界无法形成热门。

工业界的情形处于半手工作坊状态。在很多世界顶级公司中，都有相关的研究部门，数十年如一日地钻研网格生成技术；在设计一线，公司都专门培养的深谙此术的资深工程师，通过直觉洞察来手工调整软件生成的网格。这项专门的技术也通过手耳相传，世代积累，但是依然没有形成完备的理论和算法体系。工业界数十年耗费巨量资源没有解决的问题，其难度可想而知，并且肯定是在基础理论层面，而非工程层面，具有本质的障碍。

我们团队经过多年的研究和积累，逐渐认识到网格生成的基础理论不但与拓扑、黎曼几何紧密相连，而且更为关键的思想来自于共形几何和代数几何（例如全纯线丛的示性类理论），其理论基础，特别是计算理论（例如离散曲面Ricci流理论，2018）其实是近几年才发展成熟。下面，我们以曲面网格生成为例，简单介绍这方面的近期进展。

曲面三角剖分

传统的有限元方法关注在欧氏空间区域中求解椭圆型偏微分方程。例如将区域进行三角剖分，为每个胞腔构造函数基元，基函数的线性组合得到有限元的函数空间。将偏微分方程用离散或者变分法转化，在有限元函数空间中进行能量优化，从而得到代数方程。求解代数方程，得到弱解。将三角剖分加密，得到一系列弱解。如果三角剖分满足特定条件，则系列弱解收敛到光滑解。这里，三角剖分的质量至关重要，它关系到弱解序列的收敛性, 代数方程的数值稳定性等等。三角剖分的质量主要由两个指标来衡量，一个是三角形的形状，一个是三角形的尺寸。通常我们要求三角形尽量接近等边三角形，极小角有严格正下界；三角形顶点的采样密度与计算精度要求相关。经典的计算几何理论提供了成熟的算法，例如Chew算法可以保证三角剖分最小角大于,三角形大小均匀；Rupert's Delaunay Refinement算法保证最小角大于等等。

图1. 平面三角剖分的Rupert Delaunay Refinement算法。

问题的本质困难在于将平面算法推广到具有复杂拓扑和几何的曲面，这时共形几何成为解决问题的核心理论工具。更加著名的单值化定理，任何带有黎曼度量的曲面，都可以共形地映射到单位球面（用球极投影映射到平面）、欧氏平面和双曲平面上，从而将曲面网格生成转化成平面网格生成。我们先在球面、欧氏平面或者双曲平面上用经典Delaunay Refinement算法生成高质量平面网格，再将平面网格映射回到曲面上面，因为映射保角，因此曲面上的网格质量得以保证。

图2. 通过保角的黎曼映照，曲面三角剖分转换为平面三角剖分。

那么，如何计算单值化定理成为曲面网格生成的关键，同时这里有一个潜在的陷阱：离散曲面Ricci流算法的输入曲面，在计算机上输入曲面是用三角网格来表示的，这岂不是循环论证，为了生成网格我们需要输入网格？其实，问题的关键在于我们的离散Ricci流理论【2】可以保证算法对于输入曲面的三角剖分没有任何要求，无论输入网格质量如何低劣，都不会影响Ricci流算法，从而得到保角变换映射到标准曲面，标准曲面上的高质量三角剖分得到了初始曲面上的高质量三角剖分。这里离散曲面Ricci流的计算并不是传统的有限元方法在曲面上的直接推广，而是基于几何变分原理而发展出来的新颖方法，因此对于初始网格质量不敏感。由此可见，一般曲面上的非结构网格生成理论关键在于离散单值化定理，和曲面Ricci理论【2】。

图3. 用Ricci方法计算的曲面单值化。

结构化四边形网格生成

在计算机辅助设计（CAD）领域，几何曲面都是用分片有理多项式表示，即所谓的样条曲面（NURBS，T-NURBS）。样条曲面需要将曲面进行结构化四边形剖分。另一方面，结构化四边形网格剖分会提高计算精度，因此工程应用领域有大量需求。与非结构化的三角剖分相比，结构化四边形网格生成具有更大的挑战性。

多年来，学术界一直尝试建立各种理论，开发各种算法来求解四边形网格剖分。例如在计算机科学、计算数学、计算力学等领域有很多相关论文，近期比较热门的方法都是基于标架场（Frame Field）。其实，从理论角度而言，标架场依赖于曲面的微分结构和黎曼度量结构，而四边形网格剖分根本上由曲面的共形结构所决定。因此传统的标架场方法并没有揭示问题的本质，更多地是基于直观的唯像类比。

图4. 平面区域的四边形网格，奇异点构型 “牵一发而动全身”。

图4显示了平面区域的四边形剖分，每个面接近正方形，流线光滑流畅。正常的顶点拓扑度为4，与4个面相邻；红色的顶点拓扑度不等于4，被称为奇异点。经过奇异点的流线被标记成红色，被称为奇异轨道。所有奇异点的位置和拓扑度被称为奇异点的构型。我们可以想见四边形网格的整体布局取决于奇异点的构型，可谓牵一发而动全身。而奇异点之间并不彼此独立，而是相互牵制，从而具有某种整体的性质。

图5. 拓扑轮胎上不存在只有两个奇异点（3,5）的四边形网格。

通过考察更为简单的一个例子，我们可以更加强烈地感受到某种全局的障碍。众所周知，我们可以将亏格为的封闭曲面（拓扑轮胎）进行四边形网格剖分，从而不存在奇异点。一个长期悬而未决的问题是拓扑轮胎是否存在一个四边形网格剖分，只有两个奇异点，一个拓扑度为，另外一个拓扑度为？实际上，这种四边形网格并不存在，但是其证明绝非初等。如果我们假设这样的四边形网格剖分存在，并且记为，那么满足欧拉公式：点+面-边=曲面欧拉示性数，因此的存在性与拓扑不矛盾，这意味着这个问题的本质超越拓扑范畴。那么我们假设每个面都是平面欧氏正方形，如此得到曲面的一个平直度量，被称为是由四边形网格所诱导的黎曼度量。度量具有两个锥奇异点，拓扑度为的点成为一个圆锥面的顶点，锥角为，拓扑度为的点成为另一个圆锥面的顶点，锥角为。这个黎曼度量满足经典的Gauss-Bonnet定理，即总Gauss曲率等于与欧拉示性数的乘积。更进一步，根据我们上面提到的离散曲面单值化定理，我们可以证明这种带有锥奇异点的黎曼度量的存在性。因此，我们依然无法用初等黎曼几何的方法排除的存在性。

图6. 四边形网格度量诱导的和乐群。

真正揭示问题本质的观察来自度量的和乐群（holonomy）。我们在四边形网格上任选一条由面组成的回路，在第一个面中放入一个单位标架，坐标轴与四边形的两组对边平行，让后平行移动到下一个面。如此沿着回路平行移动，直至回到初始的面。这样得到的平行移动后的标架与初始标架之间相差一个旋转，转角为的整数倍。由此得到的度量曲面的和乐群是

的子群，这被称为是四边形网格的和乐条件。我们用离散曲面Ricci方法得到的平直度量虽然满足Gauss-Bonnet定理条件，但是未见得满足和乐条件，这成为问题的关键。那么，如何刻画四边形网格诱导的黎曼度量的和乐条件呢？答案也是来自共形几何。

网格诱导亚纯四次微分

给定曲面上的一个四边形网格，诱导一个自然的共形结构，即共形图册，记为，这一共形结构与黎曼度量兼容。这里兼容的意思是指，在任意的局部复坐标系上，黎曼度量满足公式：

图册可以如下建构：每个面视为复平面上的一个正方形，定义局部坐标；每条边视为复平面上的一个长方形 , 定义局部坐标；每个正常顶点，视为复平面上的一个单位开圆。这些局部坐标直接的变换都是平面上的刚体变换，具有形式

自然是双全纯的。关键是考察奇异点的邻域。假设是一个拓扑度为的奇异顶点，其邻域等距为一个圆锥面，是圆锥的顶点，锥角为。局部坐标通过指数映射将锥面展平到平面圆盘。假设是与相邻的一个面，与锥面是等距的。不妨设落在平面正方形的左下角，那么局部坐标之间的变换为：

同理，容易看出，所有局部坐标之间的变换都是双全纯的，得到了共形图册，成为一个黎曼面。

更加本质地，四边形网格不仅诱导了共形结构，并且诱导了一个亚纯四次微分，四边形网格奇异点构型等价于的极点（拓扑度小于的奇异点）和零点（拓扑度大于的奇异点）。想法非常直观，每个面上构造全纯1-形式，每条边上构造全纯微分，由局部坐标变换公式(1)，我们得到

这意味着如此构造的全纯1-形式不是全局定义的。考察奇异点邻域，由局部坐标变换公式(2)，我们得到

但是，如果我们考虑公式（3）的次幂，则

考虑公式（4）的次幂，得到

由此我们得到是亚纯四次微分，并且是全局定义的. 同时，由公式(5)我们看到拓扑度为的奇异点成为的阶极点，度为的奇异点成为的阶零点。由此，四边形网格奇异点的构型，等价于亚纯四次微分的零极点，即除子，记为

这里是在点的阶。

Abel定理

在黎曼面上，亚纯微分的除子满足著名的Abel定理，即除子满足特殊的微分、积分方程。

图7. 典范基本群基底。

假设是亏格为的封闭黎曼面，其基本群的典范基底为

黎曼面的全纯1-形式基底为

由此得到黎曼面的周期矩阵

这里

构成格点

图8. 全纯1-形式基底。

黎曼面的Jacobi簇定义为

图5中帧红色平行四边形显示了拓扑轮胎的Jacobi簇，即平环（flat torus）。Jacobi簇的点在加法下成群。Jacobi映射将黎曼面全纯地嵌入到Jacobi簇中：令，我们任选一条路径连接基点和，，

定理（Abel）黎曼面上亚纯函数的除子被Jacobi映射到Jacobi簇的零点:

我们假设黎曼面亏格为,给出Abel定理的一个简单证明。我们首先证明一个引理：

引理假设黎曼面为拓扑轮胎，是一个从球面到黎曼面的全纯映射，那么将映到上的一个点。

我们将映射提升到万有覆叠空间，为全纯映射。拓扑轮胎的万有覆叠空间与复平面共形等价，我们得到全纯映射。球面为紧集，像集有界，由刘维尔定理，全纯映射为常值映射，，这里为投影映射，因此为常值映射。

现在，我们证明Abel定理。给定一个亚纯函数，可以视为从黎曼面到球面的全纯映射，固定一点，的原像是一个除子，Jacobi映射将其映射到平环上的一点, 由此我们得到全纯映射。由引理，为常值映射。由此我们有：

这就是Abel定理。这个初等证明可以推广到任意亏格的黎曼面。

曲面四边形网格基本定理

由Abel定理，我们得到四边形网格奇异点构型的基本定理。

图9. 曲面四边形网格等价于亚纯四次微分，奇异点构型“牵一发而动全身”，表达为Abel-Jacobi条件。

定理(四边形网格奇异点构型【3】)：给定四边形网格，其诱导的亚纯四次微分为，令为任意的一个全纯1-形式，那么是一个亚纯函数，

现在我们回头用基本定理来证明轮胎曲面上不存在3-5奇异点的四边形网格剖分：假设是四边形网格，一个拓扑度为3的奇异点和一个拓扑度为的奇异点；轮胎曲面上存在一个全纯1-形式，没有零极点，除子。由上面的基本定理，我们有

这意味着点和重合，矛盾。由此，轮胎上不存在3-5奇异点的四边形网格。读者们可以想一下，轮胎曲面上是否存在一个三角剖分，每个顶点拓扑度为，但是有两个奇异点，一个拓扑度为，另一个拓扑度为？

图7. Buddha曲面转换成T-样条。

有了基本定理之后，曲面四边形网格生成的手工经验性算法可以被提升为理论严密的自动算法。通过优化奇异点的构型，使之满足基本定理的Abel-Jacobi条件(6), 然后用曲面Ricci流方法构造平直度量具有锥奇异点，从而得到亚纯四次微分，构造四边形网格，或者T样条。

小结

从现代观点来看，Abel-Jacobi定理本质是黎曼面上全纯线丛的示性类理论，亚纯微分可以看成是全纯线丛的全局截面，其除子构成了线丛的示性类。这一理论到高维复流形的自然推广就是代数几何中著名的小平邦彦定理【1】。

我们看到，工业软件的最具挑战性技术之一是曲面、实体网格生成，这项技术既有工程难度更有理论难度。工业界数十年的研发一直难以取得突破性进展，因为其理论基础需要用到非常现代的微分几何、共形几何和代数几何的核心定理，而相应的核心算法（例如离散曲面Ricci流【2】）也是近几年才被发展出来。曲面的结构化四变形网格具有强烈的几何特性，接触这一问题的所有人都立刻会产生“牵一发而动全身”的直觉。我们的工作实质上是将这一直觉用严格的代数几何的语言来表述，从而将经验性的算法提升为严密自动的算法。

目前，“数值计算引擎”非常成熟，经典的线性代数算法、数值微分方程算法和优化算法都已经成为Matlab的简单命令；但是“几何引擎”依然在襁褓之中。年轻学子们试图计算曲面两点间的测地线，曲面的调和叶状结构，或者覆叠空间，单值化度量等等，根本找不到便捷的计算软件。我们希望未来“几何引擎”能够被开发出来，微分几何、代数几何、共形几何理论方法很快在工业界推广开来，真正实现生产力的提升。

【1】Shing-Tung Yau, 'From Riemann and Kodaira to Modern Developments on Complex Manifolds', ICCM Notice, 2017.

【2】X. Gu, F. Luo, J. Sun and T. Wu, 'A Discrete Uniformization Theorem for Polyhedral Surfaces I', Journal of Differential Geometry (JDG), Volume 109, Number 2, Pages 223-256, 2018.