雪山论剑 - 工业软件发展峰会 IGA2022
横贯北美大陆的落基山脉一直延绵到加拿大阿尔伯塔省的卡尔加里。这里遍布冰川、冰原、松林、瀑布、冰湖和高山，被加拿大政府建立为国家公园。这里冬季冷峻严酷，宏伟磅礴，夏季风光旖旎，景致怡人。
Fig1. Banff 国家公园 ，Morant’s curve。（摄影 赵晗）
班芙的冰川湖举世闻名。梦莲湖折射沉积的岩粉，呈现出独特的蓝绿色，如梦如幻，宛若人间仙境。加拿大纸币的背面风景画即为梦莲湖。露易丝湖三面雪山环抱，湖水来自冰河融雪，碧绿如翡翠。这里的气候非常独特，湖边极度干燥，却又彻骨寒冷，气温在零下数十度。双手裸 露在空气中拍摄，坚持不到十秒，就已经痛彻骨髓。但是湖水清澈，景致清幽，脱离凡尘，宁静旷远。
Fig 2. 露易丝湖。（摄影 王文艳）
在原始荒野雪山之中，矗立着一座中世纪风格的古堡，班芙温泉酒店（Fairmont Banff Springs）。这家酒店以其特有的豪华与遁世，成为落基山脉的一大标志。建筑沿袭了苏格兰城堡风格，以巨大原石垒砌而成。国际等几何分析2022年会（Isogeometric Analysis 2022）在此举行。
Fig 3. Fairmont Banff Spring Hotel.
室内是中世纪的装饰，穹顶石壁，落地长窗，彩绘玻璃，烛枝吊灯，墙上是重装铠甲、骑士长矛。建筑内部回廊环绕，密道曲折，暗室嵌套，宛若迷宫。夜晚，古堡外寒风呼啸，漫天飞雪，古堡内灯光幽暗，风琴回荡。众宾客面色苍白，品味着血色红酒。特别是这里的一道招牌菜，惊悚黑暗，Alberta Tartare，即澄黄的生鸡蛋凉拌阿尔伯塔生牛肉馅，肥美多汁，鲜血迸溅，宛若德古拉的盛宴。
Fig 4. 古堡内部。
由于疫情的缘故，IGA年会停办两年，这次是疫情后的首次会议。来自北美、欧洲和亚洲学术界和工业界的专家学者们欢聚雪山之巅，切磋比试，评判研究，共同探讨工业软件和IGA发展的未来！
传统有限元方法（FEM）
计算机辅助设计（CAD）和计算机辅助工程（CAE）是工业软件的主体。CAD帮助设计者设计各种工业机械零件，用样条曲面表示各种几何曲面。CAE对设计的零件进行数值模拟仿真，即在几何实体上求解各种偏微分方程。其基本逻辑是：物理世界的规律归结为偏微分方程，工业设计产品的性能可以通过求解偏微分方程进行模拟和预测。传统CAE主要是基于有限元（FEM）方法，将几何曲面或者几何实体进行三角剖分，每个单元为单纯形。在单纯形上用重心坐标构造基函数，将函数在顶点的取值线性推广到单纯形内部，得到整体的分片线性函数。固定三角剖分，所有的线性分片函数构成一个函数空间，我们称之为有限元函数空间。将偏微分方程改写成能量变分，在有限元函数空间中进行优化，从而将无穷维变分变为有限维优化，微分方程变换成代数方程，再用线性系统求解。
由此可见，CAD领域的底层数据结构是样条曲面，例如Bezier，BSpline，TSpline，NURBS，TNURBS等等，而传统CAE领域的底层数据结构是三角网格（mesh）。这两种数据格式的转换极其繁琐和困难，通常占据整个计算流程70%以上的时间。网格生成领域的权威，德国维尔斯特拉斯研究所的资深科学家斯杭博士开发的TetGen被国际上各大CAD/CAE公司广泛采用。据斯博士介绍，网格生成的困难来源于两个方面：首先在手动CAD建模过程中，基本几何曲面的布尔运算会产生大量的裁剪样条曲面（trimmed Splines），和大量的样条碎片（feature）。一般的样条曲面是双三次多项式曲面，其交线是高次多项式曲线，无法用双三次基底函数表达，因此会产生大量的肉眼难以分辨的裂缝。自动找到这些裂缝和碎片，然后在三角剖分中将它们弥合并非易事；其次，在有限元计算中，三角剖分的质量至关重要，如果采样密度不够，计算误差会过大；如果三角形内角过小，所得刚度矩阵条件数过大，计算不稳定，即使将剖分细化，所得结果也不会收敛。如何在具有复杂拓扑和几何的曲面上生成高质量网格，一直是CAD、CAE领域的痛点。斯杭博士与笔者团队合作，用计算共形几何（Computational Conformal Geometry）来解决这个问题。一般而言，平面区域的高质量三角剖分是已经解决的问题；而共形几何可以将所有曲面保角映射到平面区域上，从而将平面上的三角剖分拉回曲面。由于映射保角，三角剖分的质量得以保持。应用笔者团队发展的离散曲面Ricci流理论，和斯杭博士的高效平面三角剖分算法，我们可以得到整体水密的曲面网格生成算法。
Fig. 5 具有复杂拓扑和几何的机械零件，由上千张裁剪的样条曲面构成。
Fig 6. 复杂机械零件用Ricci流整体共形映射到平面区域，变成一张曲面，从而得到水密的模型。
等几何分析的思想（IGA）
为了统一CAD和CAE两个领域，美国计算力学泰斗，科学院、工程院和艺术科学院的三院院士，UT Austin大学的Tom Hughes教授，在2005年左右提出了等几何分析方法（Isogeometric Analysis）。这种方法的核心是在CAE中，求解偏微分方程的时候，未知函数和测试函数所在的函数空间用样条函数空间，与输入几何曲面、几何实体的样条函数空间相一致，从而从根本上避免了网格生成问题。据Ansys资深科学家王文艳博士，Ansys的科学家李利萍博士，UCSD的赵晗博士等人的介绍，理论上，等几何分析方法具有很大的优越性：
最终实现CAD、CAE的统一，IGA使用与CAD相同的几何模型进行计算，减少前处理工作。例如省略网格生成步骤，节省70%的计算时间；
有限元方法（FEM）中生成的mesh对于输入几何实体进行逼近，具有误差；等几何分析方法（IGA）中，几何表示更加精确，误差为零；对于复杂的模型，可以精准建模；
与FEM相比，在同样精度下，IGA方法需要更少的单元计算，提高计算效率。总体而言，从每个自由度提供的精确度而言，IGA方法优于FEM。
FEM方法单元之间的连续阶比较低，IGA单元之间具有高阶连续性。而很多偏微分方程需要基函数具有高阶连续性，例如双调和方程。在这种情形，IGA成为首选。
IGA方法得到的频谱分析更加优越。
通过Bezier Extraction技术，FEM方法积累的代码可以直接用到IGA方法中，从而简化IGA的算法实现；
纽约州立大学石溪分校的陈士魁教授认为，对于剧烈形变的力学问题，传统有限元方法需要采用动态网格，计算方法相对困难，而IGA方法可以用同样的样条基底表示大形变，具有优势。例如吉林大学汽车仿真与控制国家重点实验室的陆善彬教授团队应用IGA和Ricci流方法，解决了汽车冲压成形领域的一个多年难题【5】。这里需要曲面的大形变，IGA方法更加适用。
等几何分析的发展瓶颈
虽然IGA思想先进，理论优越，但是其发展却不尽如人意，关于其内在原因，众位博士都有各自的亲身感受。
王博士、斯杭博士、Brigham Young 大学的 Kendrick Shepherd 教授认为传统的CAD方法很少考虑CAE的应用，因此构造的样条曲面无法直接应用于CAE的分析任务。例如，传统汽车行业中，汽车模型中的B-柱结构往往由上百片trimmed样条曲面拼接而成，虽然肉眼看不出瑕疵，但是物理上并非完整的一片曲面，在数值模拟中无法直接用于计算。
关于这个问题，有两个思路，一个思路是赵晗博士团队的方法，设计特殊的惩罚函数，使得曲面片接缝处零阶、1阶连续，使得接缝处网格结构错开的壳类结构可以直接应用于IGA。
另外一个思路是构造整体拼接、全局光滑的样条曲面。满足单位分解，局部可以加细，容易拓展到体的样条。这需要构造结构化四边形网格剖分。大连理工罗钟铉教授、雷娜教授团队的工作【2-4】，Hughes博士，Shepherd博士 ，笔者与Beta CAE的合作【1】就是采用这一思路。经过长期研究，我们建立了严格的数学理论，透彻地理解了这一问题的本质。这次大会，笔者作为合作者的两个演讲就是关于这一问题的解决。
王博士指出IGA的计算需要构造体样条函数，而绝大多数的CAD模型只有边界曲面的样条表示，将曲面样条拓展成内部的体样条具有根本的重要性，而这一问题是IGA发展的核心困难之一。其实，多年之前，Hughes博士就认识到曲面样条到体样条的拓展是IGA的基石。他向T-样条之父，ACM终身成就奖得主Tom Sederberg博士求助。Sederberg博士联系了笔者，因为他们认为这不是一个计算力学领域的问题，而是一个几何与拓扑问题。笔者专程飞到Austin，从Hughes博士那里学习到这一问题的背景和重要性。这一问题等价于如何将一个曲面的四边形网格剖分拓广到内部的六面体网格剖分，即网格生成领域的圣杯问题。工程界发展出了大量的经验性算法，例如将曲面变形成polycube，即长方体的并集，但是这种方法将所有的奇异点推到边界曲面，而边界的计算最为敏感，理想情形应该将奇异点深埋于实体内部。或者将实体嵌入在欧氏空间之中，在边界处微小形变成polycube。这种方法依然存在类似问题。笔者和大连理工团队发展出基于曲面叶状结构和全纯二次微分理论，构造六面体剖分的方法具有一定的普适性。但是一般曲面四边形网格是亚纯四次微分，如何将现有方法进一步推广，目前依然在探索之中。
王博士指出样条的奇异点问题。经典的样条理论是基于仿射几何不变量，而拓扑复杂曲面本身不具有仿射结构，样条曲面不可避免地具有奇异点。在奇异点，曲面的二阶连续性退化为一阶光滑性。如何构造奇异点附近具有一阶连续性的基函数，并且将这种方法拓展到体样条。这个专题的顶级学者是Ulrich Reif博士和Jorg Peters博士，Ulrich Reif博士的D-patch方法和Jorg的奇异参数化方法广为人知。Ulrich在演讲中非常推崇Jorg的方法，并且拼命阻止大家使用自己的方法，因为这种方法数学上并不完美。笔者询问Jorg，Jorg却对Ulrich的方法赞不绝口。Ulrich的演讲指出这些方法向体样条推广是会遇到本质的困难，并且告诫大家：不要被数值模拟的所谓结果欺骗，要追求理论的严密性和完备性。总体而言，样条的奇异点问题尚未严格解决。
李博士认为IGA方法对于计算断裂问题相对困难。因为每个点都由多个控制点共同支撑，因此表达非连续的断裂结构不如FEM方法灵活。
这些问题都具有本质的难度，特别是结构化曲面四边形网格剖分到结构化六面体网格剖分的拓展。Hughes博士不禁感慨：国家实验室资助世界学者数十年攻克这一问题，依然困难重重。
Fig 7. 表面T-样条拓广成体T-样条。（王文艳）
等几何分析的里程碑
Beta CAE的Lambross Rorris 博士在演讲中提到“汽车碰撞是模拟仿真的王”，（Crash is King in Simulation）. 多少年来，Hughes博士一直希望IGA方法能够真正做到汽车碰撞的模拟仿真。Kendrick Shepherd博士是Hughes博士的学生，五年读博期间跟随笔者学习计算共形几何。今年，我们终于用Ricci流方法与Abel-Jacobi条件，与Beta CAE合作，实现了汽车碰撞的分析，完成了Hughes博士的夙愿【1】。Hughes博士兴奋异常，为Kendrick颁发了银杯，因为金杯已经被Honda汽车公司于2019年夺走。
Fig 8. Crash is King in Simulation （By Lambross Rorris.）
Lambross将这一任务解释为：传统的样条曲面缺乏水密性，很多是经过裁剪的高阶多项式曲面，每一个机器部件都是由多片样条拼接而成。我们的任务是将多片拼接的传统CAD曲面转换成整体一张裁剪曲面，这样才能适合IGA的要求。而这一算法应该能够处理具有复杂拓扑和几何特征的大型几何模型，具有一定的鲁棒性和较高的效率。
Fig 9. 1996 Doge Neon Body-in-White Rebuilt Splines （Kendrick Shepherd，David Gu and Tom Hughes 【1】）
Beta CAE 提供了原始汽车模型，图9显示的是用Ricci流方法重建的样条表示，更加适用于IGA处理。
Fig 10. Crash Analysis using Conforming Splines. （Kendrick Shepherd，David Gu and Tom Hughes 【1】）
图10显示的是数值仿真的结果。我们看到碰撞后，有几个零件散落下来。Kendrick解释说，有些不重要的部分他依然用原来的样条表示，没有转换成整体连接的样条，这些曲面片具有裂缝，在模拟仿真过程中脱落下来。
Fig 11. Dodge Neon Fire Wall conversion. (K. Shepherd, D. Gu and T. Hughes 【1】)
等几何分析的几何基础
我们看到，等几何分析的基础问题本质上离不开几何和拓扑，这些理论目前还没有融入到计算力学领域。笔者的演讲主要是我们与大连理工大学团队发展出来的结构化四边形网格理论【2-4】，演讲的题目是“Structural Quadrilateral Mesh Generation Based on Abel-Jacobi Theory”。传统曲面四边形网格的构建方法主要是基于手动的经验性方法，我们的目标是为网格的奇异点构型建立微分或者积分方程，从而将算法自动化。我们的理论将曲面四边形网格与黎曼面上的亚纯四次微分建立了等价关系，用代数曲线中的Abel-Jacobi理论，建立了奇异点构型的控制方程。演讲中用到了比较多的现代几何概念，计算力学领域的学者们可能不太熟悉。台下Ulrich Reif博士询问我们的方法与传统方法的本质区别在哪里？
笔者觉得我们的理论透彻地理解了曲面四边形网格，从而指导严密算法的设计。传统的工程方法都是设计某种能量，通过优化得到能量的极值，从而得到预期结果。这种方法时而得到理想的结果，时而一无所获，这种方法本身并没有解释问题的本质，是一种唯像的经验性方法。例如，人们经常用cross field方法，就是在曲面上面覆盖光滑的cross 场。所有平面旋转，转角为90度的整数倍，构成一个4阶的对称群。位于原点的所有平面正交单位标架场的全体，构成的流形与单位圆周等价。所有正交单位标架场的集 合模掉4阶对称群，得到所有cross构成的流形，依然是单位圆周。传统的cross field方法是求从曲面到单位圆周（所有cross构成的流形）的一个光滑映射，从而得到整体光滑的cross场。这种方法理论远非严密：
曲面上的一个cross场并非是一个从曲面到单位圆周的映射，而是以单位圆周为纤维的某个纤维丛上的全局光滑截面。这句话的确切含义需要深入思考，才能理解。
上面的纤维丛并不和曲面的单位切丛等价，虽然两个丛都是同样的纤维。这里纤维丛的构造并不直观。
曲面上的cross场并不等价于四边形网格剖分。例如亏格为一的封闭曲面上不存在这样的四变形网格：只有两个奇异点，一个度为3，另一个度为5；但是同样的曲面上存在cross场，只有两个奇异点，其指标分别为正负1/4。这意味着四边形网格问题不仅仅涉及到拓扑和黎曼几何，应该也涉及到共形几何的一些深刻定理。
我们觉得，上面三点的真正理解才是彻底解决曲面样条生成的必由之路，虽然在工程领域需要时间才能够被广泛接受，但是无论哪种方法都无法绕过它们，因为这是问题的自然结构所决定的。
小结
对比于以往，这次会议吸引了很多工业界的代表：Ansys，Beta-CAE，Coreform，Sandia Lab，LS-DYNA，同时看到很多公司开始设立IGA方面的研究团队，IGA开始进入到CAE工业的主流。王博士认为IGA会在很多方面日益取代FEM。李博士认为LSDUNA是目前唯一用IGA比较成熟的商业软件，但是LSDUNA还是会保留和FEA的coupling。FEA已经发展70多年了，现在涉及的领域有板块成形，电磁，流固耦合等等。比如流体方向，即使目前为止也不是完全FEA，而是有finite volume method。可能未来各种方法 会并存，每种有它更擅长的领域而取代其他方法。
Tom Sederberg博士发明了T-样条之后成立了一家公司，由其大儿子Matt Sederberg担任CEO。很快这家公司被AutoCAD收购，同时AutoCAD垄断了T-样条的专利。这次会议上，笔者遇到了Matt，Matt目前是coreform的CEO，他非常兴奋地告诉笔者T-样条专利很快会过期，希望笔者与他们合作进行T-样条方面的算法研究。除了IGA，他们公司也非常看好拓扑优化的方向，而恰巧，笔者与陈士魁教授运用共形几何方法，将拓扑优化方法从欧氏空间推广到流形上面，并已广泛应用于共形超材料和电磁热力学结构拓扑优化设计【6-7】。陈士魁教授与笔者指导的一名希腊裔博士生【7】，毕业后成为nTopology的技术骨干。依随金属3D打印技术的成熟，拓扑优化正在如火如荼地发展。
由于疫情的原因，中国大陆的学者整体缺席，唯一的主题演讲就是笔者与大连理工团队合作的结构化四边形网格Abel-Jacobi理论。中国台湾省也有学者参加。Tom Hughes博士告诉笔者，目前全世界对于CAD、CAE领域最为重视，投资力度最大的国家就是中国。他由衷希望更多的中国学者参加明年的在法国里昂召开的IGA大会。希望到那时，IGA的根本几何拓扑问题能够得以进一步突破，IGA能够更加广泛地应用于工业界。
【1】Kendrick M. Shepherd, Xianfeng Gu and Thomas J.R. Hughes, Feature-aware reconstruction of trimmed splines using Ricci flow with metric optimization, 115555, CMAME, DOI:10.1016/j.cma.2022.115555, September 2022.
【2】Xiaopeng Zheng, Yiming Zhu, Wei Chen, Na Lei, Zhongxuan Luo, Xianfeng Gu. Quadrilateral Mesh Generation III : OptimizingSingularity Configuration Based on Abel-Jacobi Theory. CMAME 2021.
【3】Na Lei, Xiaopeng Zheng, Zhongxuan Luo, Feng Luo and Xianfeng Gu, Quadrilateral Mesh Generation II: Meoromorphic Quartic Differentials and Abel-Jacobi Condition, CMAME, 366(2020), 112980.
【4】Wei Chen, Xiaopeng Zheng, Jingyao Ke, Na Lei, Zhongxuan Luo; Xianfeng Gu, {\bf Quadrilateral Mesh Generation I : Metric Based Method}, Computer Methods in Applied Mechanics and Engineering, Volume 356, Pages 652-668, 2019.
【5】Q. Jia, X. Song, M. Ji, H. Chai, S. Lu, Isogeometric algorithm for one-step inverse forming of sheet metal. CMAME 2022.
【6】Qian Ye, Yang Guo, Shikui Chen, Na Lei, Xianfeng Gu, "Topology Optimization of Conformal Structures on Manifolds Using Extended Level Set Methods (X-LSM) and Conformal Geometry Theory", CMAME, 344 (2019): 164-185.
【7】Panagiotis Vogiatzis, Ming Ma, Shikui Chen and Xianfeng Gu “Computational Design and Additive Manufacturing of Periodic Conformal Metasurfaces by Synthesizing Topology Optimization with Conformal Mapping”, CMAME, 328 (2018): 477-497.
