用于表面多边形网格平均曲率流的域分解方法

《Computer Aided Geometric Design》：Domain decomposition for mean curvature flow of surface polygonal meshes

　　伦卡·普塔奇科娃（Lenka Ptá?ková）| 米哈尔·奥特拉塔（Michal Outrata）
捷克共和国布拉格8区索科洛夫斯卡街49/83号（Sokolovská 49/83），查尔斯大学（Charles University）数学与物理学院数值数学系（Department of Numerical Mathematics, Faculty of Mathematics and Physics）

**摘要**
我们探讨了域分解技术在处理表面网格平均曲率流（mean curvature flow）中的应用，这些网格的面是由任意简单多边形构成的。首先，我们测试了传统域分解方法在分解域是否重叠的情况下的效果。接着，我们提出了针对平均曲率流优化的罗宾耦合条件（Robin coupling conditions）的Schwarz方法。最后，我们从形状质量和纹理变形的角度分析了平滑效果。通过将初始网格分解为两个子网格，我们解决了两个较小的边界值问题，从而可以利用大量的并行计算能力来加速处理过程。

**1. 引言**
表面网格的平均曲率流（Mean Curvature Flow, MCF）是指在给定边界条件下，网格表面向所谓的“最小曲率面”演变的过程——这种曲率面的所有内点处的平均曲率均为零。MCF可以用特定的非线性偏微分方程（partial differential equation, PDE）来描述，具体内容将在第2节中详细说明。在MCF过程中，网格的平滑度是一个重要特征：随着网格的演变，高频（振荡）模式会逐渐减弱，最终只剩下低频（平滑）模式。实际上，MCF的这种平滑特性常被用于网格修复，即消除之前操作（如插值）留下的瑕疵。然而，这带来一个权衡：我们希望在平滑高频噪声的同时，保留网格的原有特征。MCF在计算机图形学中有许多应用，使其成为一个实际意义重大的研究课题。

我们针对二维离散流形（即表面网格）的MCF问题，采用了在线性偏微分方程中常用的方法体系。我们的方法基于所谓的Schwarz波形松弛方法（Schwarz waveform relaxation methods，参见Gander和Lunet 2024年第3章），这些方法结合了域分解技术（Domain Decomposition Methods, DDM，用于椭圆型PDE）和并行时间处理技术（Parallel in Time Methods, PinT，用于抛物型和双曲型PDE）。DDM/PinT方法的基本思想是将问题分解为多个子问题并并行求解，从而加快求解速度。由于这些子问题本身也是PDE，它们会因人为划分而缺少部分边界/初始条件数据。为此，我们使用初始猜测值，并通过迭代过程在子问题之间交换数据。所交换的数据类型（即耦合条件，Coupling Conditions, CC）决定了方法的类型。例如，交替Schwarz波形松弛方法交换子界面上的狄利克雷迹（Dirichlet traces）；而优化后的Schwarz波形松弛方法则交换罗宾迹（Robin traces）（具体细节见第3.4节）。

我们研究了DDM/PinT类方法在三维空间（R3）中处理多边形网格MCF过程中的应用。具体来说，对于需要平滑的表面网格，我们使用这些方法处理相应的PDE，并分析迭代方案的效率。从DDM/PinT领域的角度来看，这里的微分算子在求解过程中会发生变化（因为它受到待平滑网格的影响）。据我们所知，这是DDM/PinT领域中尚未涉及的新类型问题。从计算机图形学的角度来看，由于信息交换不完美，DDM/PinT全局近似结果可能是不连续的；只有当方法收敛到解时，才能获得更平滑的近似结果。不过，由于我们关注的是MCF过程本身而非其极限情况，因此必须调整耦合条件以确保每次迭代后的连续性和所需的平滑度。

**主要贡献**
1. 将标准子问题间耦合条件（CC）适配到非线性MCF环境中，包括为一般多边形网格定义新的法向量导数；
2. 进行数值分析（见第3节），以验证我们的方法适用于多种离散拉普拉斯算子。

**数值测试**
我们测试了多种离散拉普拉斯算子，以证明该方法不依赖于特定的拉普拉斯离散化方式。我们从多个角度评估了平滑效果，包括减少域分解引入的负面影响（如伪影）。

**结论**
尽管使用DDM/PinT方法处理MCF看似自然，但据我们所知，这一领域尚未得到充分研究。虽然DDM本身已被广泛研究，但在网格平滑方面的应用较少；同样，利用平均曲率流进行网格平滑在计算机图形学中是一个重要且研究深入的课题，但与DDM/PinT研究的联系相对有限。在介绍我们的框架之前，我们在下一节回顾了相关文献。

**1.1 相关文献**
**最小曲率面**
最小曲率面是指所有内点平均曲率为零的曲面。证明给定边界曲线下的最小曲率面的存在性被称为Plateau问题，该问题以比利时物理学家J. A. Plateau的名字命名，他通过大量实验研究肥皂泡现象。最早且极具影响力的研究之一是Pinkall和Polthier（1993年）的论文，他们提出了一个三角形网格模型，使得每个内顶点的拉普拉斯算子值为零，同时边界点保持固定或沿边界切向移动。他们使用了所谓的cotan-Laplace公式来处理离散拉普拉斯算子。

**离散平均曲率流**
离散时间下的表面网格平均曲率流传统上用于网格平滑，常用于消除不规则插值数据带来的粗糙特征。例如，我们希望去除高频噪声（即网格的高频成分），同时保留初始网格的几何特征（低频成分）。我们采用了Desbrun等人（1999年）的方法，即使用隐式欧拉方法对MCF PDE进行离散化，并用当前网格状态下的拉普拉斯算子替换原始公式中的空间非线性拉普拉斯-德拉姆算子（详见第2节）。我们将该方法扩展到任意多边形网格，并使用DDM/PinT方法进行求解。

**域分解与并行时间处理技术**
DDM/PinT方法最初作为存在性和唯一性证明的理论工具（H. A. Schwarz, E. Piccard, E. Lindel?f），后来在计算领域受到关注（P. L. Lions, M. Dryja, O. B. Widlund, P. E. Bjorstad等人）。大约100年后，这些方法在数值PDE领域得到广泛应用（参见Gander 2008年及后续参考文献）。在优化后的Schwarz波形松弛方法中，我们通过狄利克雷或罗宾耦合条件在子域边界之间交换信息。众所周知，对于强制性问题（如狄利克雷迹耦合），重叠程度越大，收敛速度越快，但计算成本也越高。当重叠不可行时，我们采用狄利克雷-诺伊曼（Dirichlet-Neumann）耦合条件以获得更好的收敛性（参见Gander和Halpern 2022年）。

**计算机图形学中的域分解技术**
DDM在计算机图形学中也有应用，例如在流体流动（English等人，2013年）和弹性模拟（Malgat等人，2015年）等任务中。然而，这些基于几何结构的方法通常需要额外的多网格或Chimera网格结构。我们的方法仅使用原始网格，并采用经典的和优化的Schwarz波形松弛方法。与Sellán等人（2019年）的边界耦合方法类似，我们的方法更为直接，属于标准Schwarz方法（未优化版本）。Sellán等人研究的是四面体网格和多个子域的情况，而我们则探索了更合适的耦合条件。Zhou等人（2024年）也研究了仅交换狄利克雷迹数据的域分解方法。

**创新点**
1. 首次在网格平滑领域测试Schwarz方法的效果，并对其可靠性进行了数值验证；
2. 首次将标准及优化Schwarz方法直接应用于一般多边形网格的MCF问题。

**总结**
我们的方法与现有计算机图形学文献的主要区别如下：
1. 我们处理的是一般多边形网格，而非仅限于三角形或四边形网格；
2. 我们的方法不依赖于多网格结构或Chimera网格；
3. 我们的耦合条件完全无重叠，某些情况下性能优于传统MCF方法（见图1.2）。

**预备知识**
设{Mt}t∈R是一组在R3中随时间演变的平滑表面集合，可能有边界也可能没有边界。若Mt满足非线性抛物线方程?2x/?t2 = H→(x)（x∈Mt, t∈R），则称Mt随时间沿平均曲率方向移动，其中H→(x)表示点x处的平均曲率向量（H→=Hn→=(κ1+κ2)n→，κ1, κ2为主曲率，H=κ1+κ2为平均曲率，n→为表面法向量）。用Δ表示Laplace–de Rham算子，我们回顾一下对于流形M?R3上点的局部坐标X，我们有(1)ΔX(x)=H→(x)，x∈M??M。我们还记得，在这种设置中，Laplace–de Rham算子与带有负号的Laplace-Beltrami算子是相同的。结合以上内容，平均曲率流旨在减小平均曲率的幅度，其表达式为(2)?X?t=?ΔX，x∈M??M，t∈R。值得注意的是，(2)与热方程相似，但关键的区别在于符号Δ的含义。在标准热方程中，定义标准Laplace算子的坐标是固定的，而在这里，我们寻找的是局部坐标——从这个意义上说，MCF可以被视为热方程的一个特定的、非线性的泛化。

2.1. 平均曲率流的时间离散化
在我们的设置中，一个可定向的表面被表示为一个多边形表面网格，或者等价地，表示为一个可定向的二维伪流形，其面可以是任何简单多边形，可能是非平面的和非凸的。因此，X对应于我们网格的顶点坐标。将(2)在时间上离散化，显式Euler方案表示为Vn+1=(I?dtΔ)Vn，其中dt表示离散时间步长，V0表示初始网格的顶点位置。Δ算子是由顶点Vn上的Laplace-de Rham算子的方便离散化得到的。我们回顾一下，显式Euler方案的稳定性是有条件的。在Desbrun等人（1999年）的研究中，作者指出稳定性准则要求dt≤min(||e||2)，其中‖e‖是边长。因此，对细网格进行有意义的平滑需要非常小的时间步长，因此必须执行大量的时间步长来通过MCF演化网格。关键的是，在每次迭代中，必须将Laplace算子应用于坐标向量，这需要集成网格元素。这在计算上是昂贵的。

2.2. 平均曲率流的隐式时间离散化
(2)的隐式时间离散化产生了一个数值稳定的方案，允许使用更大的时间步长。隐式Euler方案表示为Vn+1=Vn?dtΔVn+1，其中右侧的Δ算子应理解为依赖于所求解Vn+1的(2)中的算子。由于这个非线性问题会使求解过程在计算上具有挑战性，因此通常的做法是使用由Vn诱导的Δ。这样，一步仍然需要求解一组代数方程，即(3)(I+dtΔ)Vn+1=Vn，但由于Δ是由Vn诱导的，这个系统现在是线性的。

2.2. 离散拉普拉斯算子
在过去的十年中，出现了许多不同的Laplace-de Rham和Laplace-Beltrami算子的离散化方法，在MCF的背景下，这些方法的性能差异很大。计算机图形学中网格平滑的一个重要方面是保持网格的纹理——显然，任何程度的平滑都会导致纹理失真，但只要失真不明显，我们通常对结果感到满意，见图1。为了证明我们对Schwarz方法的改编不是针对特定的拉普拉斯算子离散化，我们使用了三种不同的拉普拉斯算子离散化方法；所有这些方法都是为一般多边形网格定义的，我们接下来以矩阵形式展示它们。我们首先定义一组构建矩阵，即其乘积导致离散拉普拉斯算子的矩阵。这些构建矩阵在第3.2节中也用于计算法线导数。设M是一个具有顶点V、定向半边E和定向面F的网格。我们还用?,?表示关联关系，即f?v表示v是面f的顶点，我们使用关联数[fi:ej]，即[fi:ej]=1表示ej是与fi具有相同方向的半边。由于我们使用的是半边结构，我们有[fi:ej]=1或[fi:ej]=0。我们用d0∈R|E|×|V|和d1∈R|F|×|E|表示在0-形式和1-形式上的离散外导数矩阵，分别为(4)d0[i,j]=[ei:vj]={1如果vj是ei的端点，-1如果vj是ei的起点，0否则，(5)d1[i,j]=[fi:ej]={1如果ej?fi且ej的方向与fi相同，0否则。进一步，矩阵fv∈R|F|×|V|表示面和顶点之间的关联关系，表示为(6)fv[i,j]={1如果vj?fi且fi是pi-边形，0否则。矩阵R∈R|E|×|E|，R=diags(Rf0,…,Rf|F|?1)，表示边的邻接关系，它是根据每个pi-边形fi给出的块对角矩阵，表示为(7)Rf=∑a=1?p?12?(12?ap)Ra，Ra[k,j]={1如果ej是f的第(k+a)个半边，[f:ej]=1，-1如果ej是f的第(k?a)个半边，[f:ej]=1，0否则。矩阵A∈R|E|×|E|则编码一条边是否有两个半边（它是内部边）或只有一个半边（它是边界边）。它取与同一条边相邻的定向半边值的线性组合：(8)A[i,j]={1如果i=j且ej是边界边，1/2如果i=j且ej是内部边，-1/2如果ej是与ei方向相反的同一条边的半边，0否则。注意，所有上述矩阵都是与度量无关的，因为它们只捕获网格元素的连接和关联，因此我们称它们为拓扑矩阵；它们的构建通常占相关实验和计算总运行时间的大部分。我们还定义了与度量相关的矩阵WV∈R|V|×|V|、W1∈R|E|×|E|和向量wE∈R|E|×1，表示为(9)WV[i,i]=1∑fk?vi|fk|pk，(10)W1[i,j]={?ei,ej?|f|如果?f使得[f:ei]=1且[f:ej]=1，0否则，(11)wE=(1‖e0‖2,…,1‖e|E|?1‖)?，其中(e0,…,e|E|?1)?=d0V。由于我们处理的网格的面是p-边形，可能是非平面和非凸的，我们计算面f的面积为其向量面积的模：|f|=1/2‖∑i=0p?1vi×v(i+1)mod p‖2，对于f=(v0,…,vp?1)。面f的向量面积的模|f|是f在R3中所有正交投影上的最大面积。向量面积仅依赖于边界曲线，而不依赖于支撑表面（Alexa和Wardetzky（2011年，引理1））。

Fujiwara（1995年）的尺度依赖性伞形算子由以下公式给出Δ(vi)=1/m(vi)∑vj～viivi?vj‖vi?vj‖2，m(vi)=∑vj～vi‖vi?vj‖2，其中∑vj～vi表示通过边连接的顶点vi,vj的所有对的求和。我们通过将其乘以vi的价数来调整拉普拉斯算子，使其与Gander和Halpern（2022年，第2节）在均匀网格上的标准二阶五点有限差分离散化相匹配，只是符号不同。因此，我们得到了调整后的Fujiwara拉普拉斯算子的以下矩阵形式(12)ΔF=1/2|d0?|?(1,…,1)?|d0?|?wEd0?diag(wE)d0，其中绝对值|d0?|是指逐元素。Alexa和Wardetzky（2011年）的cotan拉普拉斯算子和多边形拉普拉斯算子。原始的cotan拉普拉斯算子在Pinkall和Polthier（1993年，方程（6））中定义在三角形网格上，不仅考虑了相邻边的长度，还考虑了与vi相邻的角度。最近，在Alexa和Wardetzky（2011年）中，作者提出了一类在一般多边形网格上的离散拉普拉斯算子，推广了cotan拉普拉斯算子——见其中的定理2。我们关注他们的几何拉普拉斯算子，这里表示为ΔA，它可以表示为矩阵形式(13)ΔA=WVd0?RW1R?d0。Ptá?ková和Velho（2021年）的多边形拉普拉斯算子是通过离散外微分导数定义的，对应于离散Laplace–de Rham算子，也称为Hodge拉普拉斯算子，因为它使用他们的离散Hodge星算子定义。我们表示他们的拉普拉斯算子为ΔP；其表示为(14)ΔP=?WVfv?d1AW1R?d0。注意，所有上述拉普拉斯算子都对应于Laplace–de Rham算子，因此它们是Laplace–Beltrami算子的负值，后者对于函数g定义为Δg=?div(?g)。

2.3. 初始和边界条件
我们通过陈述伴随我们离散化MCF（3）的离散初始和边界条件来完成问题设置。初始条件V0是我们想要平滑的网格的顶点位置矩阵。此外，我们假设我们有（用户给出的）将顶点划分为在平滑过程中不能移动的边界顶点β，以及在平滑过程中可以自由移动的内部顶点ι，并且对于任何n≥0，我们按Vn=((Vβn)T,(Vιn)T)T的顺序排列Vn。然后，假设表面网格有边界，即β非空，我们从(3)中消除Vβn，而不是(I0dtΔβI+dtΔι)(VβnVιn+1)=(VβnVιn)，我们得到较小的系统(15)(I+dtΔι)Vιn+1=Vιn?dtΔβVβn，其中选择的离散拉普拉斯算子已根据索引(β,ι)分解为子矩阵。接下来，我们省略下标Δβ,Δι以避免符号混乱。我们以一个具有纹理的一般多边形网格的MCF示例来结束这一节。

示例1
在图2中，我们展示了MCF对由一般多边形形成的网格的影响。映射在表面上的纹理是原始网格的高度图。颜色图突出显示了ΔF在H→中引入了切向分量，因此纹理在过程中发生了漂移。虽然使用ΔA的流动会出现不希望的伪影，但ΔP的表现很好。另一方面，在四边形网格上，使用Fujiwara（1995年）和Alexa和Wardetzky（2011年）的拉普拉斯算子得到的平滑效果通常比使用Ptá?ková和Velho（2021年）的拉普拉斯算子更好。这也是我们测试我们的DDM与Robin CC结合这三种不同拉普拉斯算子离散化的原因之一——它们各有优缺点，根据具体意图，我们可能会选择其中之一。

3. 带有域分解的平均曲率流
我们有一个想要平滑的网格M，首先将其划分为两个相连的子网格MA和MB。这两个子网格可以重叠，或者只共享一个公共接口，见图3。然后我们计算通过适当的CC耦合的MA和MB的MCF：在隐式Euler的每一步中，我们根据MB（和MA）中的顶点为MA（和MB）指定边界条件ΓA（和ΓB）。我们期望这两个子网格保持具有边界的二维伪流形的拓扑结构，并测试了三种不同类型的CC：1. Dirichlet CC，产生类似于经典Schwarz波形松弛方法的方法；2. Robin CC，产生类似于优化后的Schwarz波形松弛方法的方法；3. 混合CC，产生类似于经典和优化后的Schwarz波形松弛方法混合的方法。我们想强调的是，所有这些方法都与启发它们的DDM/PinT方法有显著不同。首先，我们回顾一下，在MCF过程中，离散拉普拉斯算子以及法线和切线导数的模板会发生变化，这与Schwarz波形松弛（以及一般的DDM/PinT方法）的标准设置不同。其次，Schwarz波形松弛方法用于独立解决每个空间子域在整个感兴趣的时间间隔内的初始值问题，然后通过CC交换接口数据并再次求解，开始一个迭代过程（在我们的情况下这在计算上将是极其昂贵的）。我们的提议是在隐式Euler方法的每一步之后强制耦合，以便MA的MCF尽快将接口变化传播到MB的MCF。

3.1. Dirichlet耦合条件
考虑到标准热方程，使用Dirichlet CC的经典Schwarz波形松弛方法只有在有重叠的情况下才会收敛，并且收敛速度比使用Robin CC的优化版本慢，见Bennequin等人（2016年），Gander和Lunet（2024年），假设重叠是固定的；这与Gander（2006年）对Poisson问题的经典Schwarz方法的结果类似。虽然实验上这些性质在我们的设置中仍然成立，但在某些情况下使用Dirichlet CC会很有用。例如，如果网格是一个非常不规则的多边形网格，并且接口没有放置在高频模式区域，见图4中的示例。下面给出了经典Schwarz波形松弛方法的MCF改编的算法描述。用VA,B表示子网格MA,B的顶点，用ΓA,B表示MA,B的接口顶点，见图3。具有重叠的交替MCF-Schwarz方法如下：给定顶点VA,Bn，我们首先通过MCF演化MA一个时间步长，将其顶点沿ΓA固定为VBn中的顶点，即解决(16){(I+dtΔ)VAk+1=VAk,VA∈MA??MA,VAk+1=VAk,VA∈?MA?ΓA,VAk+1=VBk,VA∈ΓA,VB∈MB∩ΓA，见图3，然后通过MCF演化MB一个时间步长，将其顶点沿ΓB固定为(16)中计算的顶点，即VAn+1，解决(17){(I+dtΔ)VBk+1=VBk,VB??MB,VBk+1=VBk,VB∈?MB?ΓB,VBk+1=VAk+1,VB∈ΓB,VA∈MA∩ΓB。解决(16)和随后的(17)对应于该方法的一个步骤，交替流动MA和MB。在图4中，我们展示了具有和不具有分解的一般多边形网格的MCF。我们选择了一个设置，在该设置中无法简单地应用Robin CC，因为不清楚如何定义接口上的法线和切线导数。然而，由于接口位于一个相对平滑的区域，即我们不期望沿接口发生大量平滑，因此其顶点在MCF过程中不会移动太多，交替的MCF-Schwarz方法表现良好。然而，当我们使用这种方法处理一个界面应该发生显著移动的网格时，我们会得到不希望出现的伪影，如图8所示。在那里，界面顶点的移动速度比内部顶点慢得多，导致由于该方法而产生的可见伪影。3.2. Robin耦合条件如第3.1节所述，改进Dirichlet CC的标准方法是考虑沿界面的（高阶）导数迹。这样做的动机有两个方面。首先，我们可以探讨什么样的CC是最优的，以便Schwarz方法的解与不进行网格分解时的解一致。在线性偏微分方程的情况下，这一点已经被广泛研究过，对于热方程这个典型例子，我们得到了Dirichlet到Neumann算子。Robin CC是其在频域中的自然零阶近似，参见Gander和Lunet（2024年，第3.4节）。其次，我们可以探讨什么是改进Dirichlet CC的最简单方法，该方法在计算上不昂贵并且具有几何解释。答案是规定沿界面的Dirichlet迹和Neumann迹的线性组合——即所谓的Robin CC——而不仅仅是Dirichlet迹。如果线性组合主要由Dirichlet迹主导，那么我们接近第3.1节中的方法，虽然保证了连续性，但可能会有伪影；如果线性组合主要由Neumann迹主导，那么我们消除了这些伪影，但我们不再保证连续性。强制实现迹的适当平衡，则会得到比任一极端情况都更好的CC。此外，离散模板保持局部性，并且与Dirichlet CC相比，不会显著影响计算效率。一般来说，我们还可以在线性组合中包括切向导数，从而得到完整的Ventcell CC，参见Gander和Halpern（2022年，方程（1））。然而，如果界面高度弯曲，那么使用切向导数似乎不太明智，因此我们在这里不进一步考虑这种推广。稍微滥用适当的顺序，我们首先陈述条件，并在第3.4节中定义离散法向导数。让?n表示离散法向导数（见定义1）。如第3.1节所述，我们用VA,B表示MA,B的顶点，但这里我们考虑网格M的非重叠划分，见图3，它们有一个共同的界面Γ≡ΓA,B。优化后的MCF-Schwarz方法结合Robin CC对应于对MCF的Robin CC的直接适应，其步骤如下：给定顶点VA,Bn，我们首先通过MCF使MA演化一个时间步长，使其沿ΓA的顶点满足由VBn中的顶点给出的Robin/Ventcell CC，即解决(18){(I+dtΔ)VAk+1=VAk,MA??MA,VAk+1=VAk,VA∈?MA?ΓA,(I+dt?p?n))VAk+1=(I+dt?p?n)VBk,VA∈ΓA,VB∈ΓB，然后通过MCF使MB演化一个时间步长，使其沿ΓB的顶点满足由(18)中计算的顶点给出的Robin/Ventcell CC，即通过VAn+1解决(19){(I+dtΔ)VBk+1=VBk,MB??MB,VBk+1=VBk,VB∈?MB?ΓB,(I+dt?p?n)VBk+1=(I+dt?p?n)VAk+1,VA∈ΓA,VB∈ΓB。解决(18)然后是(19)对应于该方法的一个步骤。通常，这种方法会导致网格不连续，这使得其在网格平滑中的实际使用受到质疑。接下来，我们引入对((18), (19))的适应，以便网格在整个MCF过程中保持连续。3.3. 混合耦合条件我们调整了优化的MCF-Schwarz方法，以便在每次迭代之后，即在通过MCF演化MA和MB之后，重新组装的网格M=MA∪MB?ΓB是连续的。我们通过结合优化的MCF-Schwarz方法和交替的MCF-Schwarz方法来实现这一点，即我们使用第3.2节中的相同MCF计算MA，使用Robin CC，但使用第3.1节中的Dirichlet CC来计算MB，从而确保MA∪MB?ΓB的连续性。用?n表示离散法向导数（见定义1），我们再次用VA,B表示MA,B的顶点，用VB表示MB的顶点，并假设这些形成了网格M的非重叠划分，见图3，它们有一个共同的界面Γ≡ΓA,B。结合Robin CC的混合优化MCF-Schwarz方法如下：给定顶点VA,Bn，我们首先通过MCF使MA演化一个时间步长，使其沿ΓA的顶点满足由VBn中的顶点给出的Robin CC，即解决(20){(I+dtΔ)VAk+1=VAk,MA??MA,VAk+1=VAk,VA∈?MA?ΓA,(I+dt?p?n)VAk+1=(I+dt?p?n)VBk,VA∈ΓA,VB∈ΓB，然后通过MCF使MB演化一个时间步长，固定其沿ΓB的顶点以对应于(20)中计算的顶点，即通过VAn+1解决(21){(I+dtΔ)VBk+1=VBk,VB∈MB??MB,VBk+1=VBk,VB∈?MB?ΓB,VBk+1=VAk+1,VA∈ΓA,VB∈ΓB。解决(20)然后是(21)对应于该方法的一个步骤。定义了这些方法后，我们现在关注符号?n的精确含义，并为一般多边形网格定义了沿界面的法向导数。3.4. 离散法向和切向导数离散法向导数是在网格M上沿着“垂直”于边界Γ的方向的方向导数的离散版本。如图5所示，在顶点vi处垂直于边界曲线的方向甚至没有连续方向导数的良好定义。因此，我们在以下定义中提出了它的近似。从几何角度来看，离散法向导数是垂直于界面曲线的向外指向的切向向量的近似。定义1设vi?1, vi, vi+1是给定网格M边界Γ上的三个连续顶点。用c(fj)表示与vi相邻的面fj的质心，即c(fj)=1/pj∑vk?fjvk，对于一个pj边形ff。进一步，用t→=vi+1?vi?1‖vi+1?vi?1‖2表示。顶点vi处对Γ的离散法向导数?n为：(22)?n(vi)=N→??N→,t→?/t→‖N→??N→,t→?/t→‖2，N→=∑fj?vi(vi?c(fj))。在法向导数的定义中，我们首先总结了起点是与顶点vi相邻的面的质心且终点是顶点vi的向量。我们用N→表示这个向量，见图5。然后我们将N→投影到法向量为(vi+1?vi?1)的平面上，然后对其进行归一化；这是必要的，因为我们也在处理各向异性网格。例如，在图3中，界面上的中心顶点在左侧子网格中有两个相邻的面，在右侧子网格中有三个相邻的面。在平面均匀网格上，上述定义将对应于Gander和Halpern（2022年，第2节）中给出的标准二阶五点有限差分离散化的法向导数。在我们的实现中，我们可以使用方程（6）中的矩阵fv轻松计算面的质心：(c(f0)?c(f|F|?1))=(diag((1,…,1)?fv?)??fv??fv)(V0?V|V|?1,)其中天花板函数???是逐元素应用的。再次设vi?1, vi, vi+1是界面Γ上的三个连续顶点。顶点vi处沿界面Γ的离散切向导数?t，即方向导数为：(23)?t(vi)=2‖vi?vi?1‖+‖vi?vi+1‖/(vi?vi?1‖vi?vi?1‖+vi?vi+1‖)，其中范数是标准的欧几里得范数。同样，在均匀网格上，切向导数将对应于Gander和Halpern（2022年，第2节）中给出的标准有限差分离散化。3.5. Robin参数p的选择我们进行了大量实验，使用了不同的时间步长、三种不同的离散拉普拉斯算子，在许多不同的网格上，数值近似参数p的值，从而获得了结合Robin CC的混合优化Schwarz方法的最佳性能。我们在下面的第4节中展示了使用该方法的经验，即对于热方程的优化Schwarz波形松弛，渐近最优的参数p?选择Robin界面条件可以接近最佳性能。因此我们使用p?的封闭形式公式，见Gander（2006年，方程（4.15）：(24)p?=((kmin2+η)(kmax2+η))1/4，其中η=1，kmin是网格可以表示的最低频率，kmax是给定网格可以表示的最高频率。根据Gander（2006年，定理4.2的证明）的发现，并适应我们的情况，我们设置kmax=π‖e‖，kmin=π?，其中‖e‖是输入网格的最小边长，?是初始子网格的宽度（Γ的顶点与对面边界的顶点之间的近似距离）。总的来说，方程（24）表示为(25)p?=(((πl)2+1)((π‖e‖)2+1))1/4。尽管我们在这里没有详细考虑Ventcell CC，但我们建议感兴趣的读者参考Gander（2006年）；Bennequin等人（2016年）以提取提到的Ventcell CC的适当参数p,q——如果界面Γ上的曲线相当直，这些参数可能会很有用。例如，在网格状网格上，如果我们使用ΔF进行MCF，我们可以使用切向导数而不会引入任何伪影，如图10所示，因为Fujiwara（1995）的拉普拉斯算子本身对界面曲线有拉直效果。注释1使用结合Robin CC的混合优化Schwarz方法对MA,B进行MCF（通过解决（20），（21）），我们得到ΓA=ΓB。然而，重新组装的网格M=MA∪MB?ΓB在接近界面Γ的地方可能会变得不光滑。我们通过在每次完整迭代后进行校正来轻松解决这个问题。我们通过替换顶点vi∈ΓA为其邻域顶点的质心来进行非常简单的校正，即重新定义它们为(26)vi=1/n∑vj～vivj，其中n是vi的价数。下载：下载高分辨率图像（567KB）下载：下载全尺寸图像图1. 通过平均曲率流（MCF）对表面多边形网格进行平滑处理（最左边），其上半部分的细节水平是下半部分的两倍。我们将网格分解为两个子网格（分隔上下两部分）。然后应用五次MCF和我们的适应Robin CC，时间步长dt=0.05。结果网格分别对应于Alexa和Wardetzky（2011年）（第二张图像）、Ptá?ková和Velho（2021年）（第三张图像）以及Fujiwara（1995年）（第四张图像）的拉普拉斯算子。最后第五张图像是未经分解的Fujiwara（1995）拉普拉斯算子的MCF结果网格。我们可以看到，尽管子网格的细节水平不同且界面是一条曲线，但我们适应的Robin CC工作得很好，结果网格也非常平滑。在Fujiwara（1995）的拉普拉斯算子的情况下，我们的域分解方法实际上防止了由于细节水平不同而导致的顶点切向移动，这可以在第五张图像中观察到明显的纹理变形。下载：下载高分辨率图像（323KB）下载：下载全尺寸图像图2. 使用平均曲率流（15）对多边形网格上的凸起进行平滑处理一次，时间步长dt=0.01，使用Fujiwara（1995年）（中间左）、Alexa和Wardetzky（2011年）（中间右）以及Ptá?ková和Velho（2021年）（最右边）的拉普拉斯算子。在一般多边形网格上，Ptá?ková和Velho（2021年）的拉普拉斯算子通常会得到更好的结果，因为Alexa和Wardetzky（2011）的拉普拉斯算子很快会产生伪影，而Fujiwara（1995）的拉普拉斯算子会导致切向移动，从而导致纹理变形。下载：下载高分辨率图像（203KB）下载：下载全尺寸图像图3. 将网格M分解为两个子网格MA和MB，一个面宽度的重叠（左）和没有重叠（右）。每个网格的顶点分别存储，因此属于重叠部分（或仅属于界面的）的顶点存储在两个副本中，一个属于MA，另一个属于MB。下载：下载高分辨率图像（481KB）下载：下载全尺寸图像图4. 对一般多边形网格进行MCF处理（最左边），该网格已分解为两个有重叠的子网格（中间左）。我们应用两次MCF，时间步长dt=0.05，并使用Ptá?ková和Velho（2021）的拉普拉斯算子，首先对整个网格（中间右），然后使用方法（（16），（17）（最右边）。下载：下载高分辨率图像（181KB）下载：下载全尺寸图像图5. 顶点vi处界面ΓA的离散法向导数?n（红色）是一个向量?n(vi)，用洋红色表示。它对应于向量N→在法向量为(vi+1?vi?1)（绿色）的平面上的归一化正交投影。向量N→（青色）是向量（vi?c(fj?1)）和（vi?c(fj)）（蓝色）的和，其中点c(fj?1)，c(fj)是与顶点vi相邻的面的质心。（关于此图例中颜色的解释，请参阅本文的网页版本。）下载：下载高分辨率图像（316KB）下载：下载全尺寸图像图6. 我们适应的MCF-Schwarz方法在主要为六边形的网格上（|V|=6118，min‖e‖≈0.005）。原始网格（最左边）自动分为两个大致相等的子网格；dt=0.01。然后跟踪平滑参考网格与我们适应的MCF-Schwarz方法（（20），（21）平滑后的网格之间的相对距离，对于ΔF，ΔA和ΔP。y轴表示迭代次数，x轴显示相对距离（有关图例的解释，请参见第4.2节）。除了ΔA之外，相应网格的距离都很小。然而，对于ΔA，即使不进行分解，从第5次迭代开始，网格也会变成一个不光滑的尖锐网格。下载：下载高分辨率图像（360KB）下载：下载全尺寸图像图7. 我们改进的MCF-Schwarz方法在三角形网格上的应用（3721个顶点，min‖e‖≈0.03）。网格自动分为两个大致相等的子网格；时间步长为dt=0.05。然后观察平滑参考网格与我们方法平滑后的网格之间的相对距离图（（20），（21）对于ΔF、ΔA和ΔP。y轴表示迭代次数，x轴表示相对距离（有关解释，请参见第4.2节）。我们想补充的是，对于ΔP，即使不进行分解，从第3次迭代开始，流动也会导致网格变得不光滑且尖锐。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像图8. 对于具有10k个顶点的四边形网格，应用了域分解的MCF（第一张图）。我们进行了5次迭代，时间步长为dt=0.05，并展示了使用DDM和我们改进的Robin CC平滑后的网格（第二张图）以及ΔF（第三张图）和ΔP（第四张图）。我们还展示了使用标准Schwarz交替方法和平面大重叠平滑后的网格（第五张图）。标准Schwarz交替方法引入了显著的伪影，即界面上的顶点移动到（最终）最小表面的速度比内部顶点慢得多。下载：下载高分辨率图像（1MB）下载：下载全尺寸图像图9. 对于具有10k个顶点的不规则各向异性三角形网格，应用了域分解的MCF（第一张和第二张图）。我们应用了5次我们改进的MCF-Schwarz方法，时间步长为dt=0.05，并以伪彩色显示了结果网格与未分解平滑网格之间的H_max距离（第四张图），其中H_max = 0.12500244（深红色）。（关于此图例中颜色的解释，请参阅本文的网页版本。）下载：下载高分辨率图像（2MB）下载：下载全尺寸图像图10. 通过平均曲率流对网格进行平滑处理，无论是否进行域分解；时间步长为dt=0.05。左上角是原始网格及其高度图的纹理；底层表面与图2相同。然后从右到左依次展示使用Fujiwara（1995）、Alexa和Wardetzky（2011）以及Ptá?ková和Velho（2021）的Laplacians进行1次MCF处理后的网格。在底部行中，我们在首先将原始网格分解为两个子网格后进行平滑处理；其中一部分是最左边的图像。我们使用与上一行相同的Laplacians顺序。当我们在界面上选择适当的条件时，很难区分相应的网格。4. 实现与评估4.1. 实现我们已经使用Python编程语言实现了Schwarz交替方法（（16），（17）以及带有Robin CC的混合优化Schwarz方法（（20），（21）。代码可在GitHub仓库中找到：https://github.com/lenka-ptackova/ddm-mcf。在我们带有Robin CC的混合优化Schwarz方法中，我们在两个实例上使用了并行性：（1）用于创建编码子网格拓扑的稀疏矩阵（（4），（5），（6），（7），（8），这些矩阵在整个MCF过程中不会改变；（2）用于构建和更新对应于子网格MA和MB的离散Laplacians。我们进行了大量测试（请参见上述GitHub仓库），并提供了一个总结我们整体经验的例子。scipy代码在配备Intel Core i7-1165G7 2.80 GHz处理器和16 GiB RAM的笔记本电脑上运行；并行计算是在CPU上使用object multiprocessing.Pool完成的，每个子网格使用一个核心。在表1中，我们报告了使用Alexa和Wardetzky（2011）的Laplacians对两种不同网格进行两次不同细化处理的时间和加速比。表1. 未分解和分解情况下MCF子例程的平均时间以及相对加速比。Laplacians使用的是Alexa和Wardetzky（2011）的，但其他结果通常也是可比的。网格（|V|）MCFMCF DDM加速比拓扑Laplacianspsolve拓扑Laplacianspsolve拓扑Laplacianspsolve图8（10k）123.12秒1.293秒0.060秒32.65秒0.845秒0.069秒3.77倍1.53倍0.88倍图8（28k）952.75秒3.473秒0.224秒253.16秒3.161秒0.225秒3.76倍1.1倍0.99倍图9（10k）278.01秒1.852秒0.056秒74.42秒1.165秒0.055秒3.74倍1.6倍1.03倍图9（28k）2100.2秒5.038秒0.214秒58.858秒3.977秒0.226秒3.61倍1.27倍0.94倍我们看到，由于并行性，即使是中等大小的网格也能获得显著的总体加速。在我们看来，这证明了概念的可行性，明确鼓励进一步的实验研究，特别是在多个子网格的方向上，可能会带来巨大的加速，而不会明显影响平滑网格的质量。我们还发现，构建拓扑矩阵是迄今为止最耗时的任务，其次是更新Laplacians。尽管我们按顺序解决了（20），（21），但没有观察到明显的减速。有关并行化和标准DDM/PinT方法加速比的更多参考文献和分析，请参阅Dolean等人（2015）；Gander和Zhang（2019）；Gander和Lunet（2024）。4.2. 几何评估我们通过比较分解后流动并重构的网格M~k与未分解使用MCF流动的参考网格Mk来评估DDM在MCF中引入的误差。我们通过视觉方式以及通过测量相应顶点之间的距离和使用PyMeshLab库中的hausdorff_distance方法测量网格M~k和Mk的伪Hausdorff距离来进行客观比较。hausdorff_distance方法返回三个量：H_RMS——参考网格顶点到候选网格上最近点（不一定是顶点）的距离的均方根；H_max——参考网格的每个顶点与候选网格上最近点之间的距离的最大值；H_mean——这些距离的平均值。我们同样表示参考网格和候选网格相应顶点之间的相对距离：V_RMS——相应顶点距离的均方根；V_max——最大距离；V_mean——距离的平均值。在图6和图7的图表中，我们绘制了相对距离，即将计算出的距离除以参考网格Mk的平均边长。在图9中，我们展示了所提出的方法也适用于各向异性不规则网格；我们根据顶点与我们方法平滑后的网格之间的伪Hausdorff相对距离对网格的顶点进行了着色。5. 结论我们研究了几种耦合条件以及三种不同的Laplacian离散化方法，因为每种方法都有其优缺点。因此，根据网格的不同，用户可以选择其中之一。通过展示所提出的域分解方法与这三种离散化方法都能很好地配合使用，我们证明了所提出的框架具有广泛的适用性。具体来说，Alexa和Wardetzky（2011）以及Ptá?ková和Velho（2021）的Laplacians产生的平均曲率向量的切向分量可以忽略不计。因此，使用这些Laplacians的MCF在移动过程中不会变形映射在网格上的纹理。Fujiwara（1995）的Laplacians导致的顶点切向移动可以在图2和图8中观察到。这一观察结果也适用于带有域分解的MCF。另一方面，Fujiwara（1995）的Laplacian在规则各向同性网格上最为稳定。在第3.1节中，我们介绍了适应于MCF任务的Schwarz交替方法。我们重申，如果子网格之间的界面位于平均曲率较低的区域，或者预计在流动过程中不会显著移动的区域（例如沿着“直线”，如图4所示），这种方法是有用的。另一方面，如果界面穿过一个峰值，如图8所示，我们会观察到界面顶点的流动相对于子网格的内部顶点有显著延迟。在一般的多边形网格上，由于难以选择合适的子域界面（如图4所示），定义法线导数可能会很困难。然而，在界面上满足解释条件的情况下，可以使用经典的交替Schwarz方法并获得良好的结果。定义1中给出的Robin CC和（25）中给出的p?保证了在应用未分解的MCF时，演变表面网格的总体形状将模仿相应网格族的形状。界面不必是直线，但它必须允许定义良好的法线导数。此外，我们还展示了我们的改进方法允许在分解的情况下进行MCF，而不会引入如纹理变形之类的伪影。更重要的是，当初始网格在某些区域具有不同的细节级别时，我们的改进方法也是可靠的，如图1所示。我们还证明了使用改进的MCF-Schwarz方法由于计算中最昂贵部分的并行化，可以获得显著的加速（超过3.5倍）。总之，我们认为我们的阐述为适应和（数值上）探索使用DDM/PinT方法进行网格平滑提供了概念证明。主要是因为在并行计算时代，它们可以保证方法的加速，同时不会显著影响平滑网格的质量，正如第4节所展示的。CRediT作者贡献声明Lenka Ptá?ková：概念化、形式分析、研究、方法论、项目管理、软件、验证、可视化、撰写——原始草稿、撰写——审阅与编辑。Michal Outrata：概念化、形式分析、资金获取、研究、方法论、资源、撰写——原始草稿、撰写——审阅与编辑。