ParFEMWARP：一种基于并行有限元的网格变形算法——算法设计、分析与实现

《Computer-Aided Design》：ParFEMWARP, The parallel finite element-based mesh warping algorithm: Algorithm design, analysis, and implementation

　　•ParFEMWARP 是基于有限元的网格变形算法的并行化版本。
•我们的实现可以在共享内存/分布式内存架构上移植。
•我们通过分离物理/运动逻辑的抽象来实现自定义变形。
•我们为方法的每个阶段提供了并行运行时间和空间复杂度分析。
•我们在全局刚度矩阵生成步骤中展示了理想的加速效果，在共享和分布式内存架构上整体方法实现了亚线性加速。

**引言**
网格变形对于建模动态问题非常重要，例如生物力学、空气动力学、图像处理和心脏病学应用中的问题。具体用例包括在疾病进展和治疗过程中跟踪大脑运动[1]、飞机机翼的表面变形[2]、图像变形[3]、将多张图像拼接成一张[4]、跟踪心律[5]以及变形参考器官以创建患者特定模型[6]。尽管对来自实际应用的大网格进行变形在计算上非常昂贵，但它可以从某种形式的并行化中受益，无论是多核、节点、GPU等。因此，我们提出了由 Shontz 和 Vavasis 在 [7]、[8] 中提出的基于有限元的网格变形算法（FEMWARP）的并行化版本。该方法首先使用有限元方法（FEM）解决泊松边值问题（BVP），然后根据初始网格和边界变形解决具有多个右侧的线性系统。

几种与网格变形和 FEM 相关的并行算法在各种架构上显示出很强的扩展性。Sastray 等人开发了一种基于 OpenMP 的共享内存并行化版本 FEMWARP 用于变形四面体网格[9]。Secco 等人开发了一种基于 MPI 的分布式内存并行化版本逆距离（IDWARP）方法用于变形结构化六面体网格[10]。此外，Gerhold 和 Neumann 开发了一种基于域分解的分布式内存、MPI 并行化版本的非结构化网格变形方法[11]。Panitanarak 和 Shontz 开发了一种基于MPI的分布式内存并行化版本的对数障碍变形算法用于变形四面体网格，但每个节点使用一个核心以最大化每个等级的可用内存，从而实现强大的可扩展性[5]。Zhao 等人开发了一种基于径向基函数（RBF）的混合 MPI/OpenMP 并行化网格变形算法，并获得了强大的可扩展性[12]。Luke 等人开发了一种基于逆距离加权函数的并行网格变形算法，用于非结构化多面体网格[13]。他们的方案产生的网格质量与 RBF 方法类似，但在保持计算流体动力学应用中粘性网格的边界层正交性方面表现更好。Krysl 使用共享内存和 pthreads 实现了稀疏 FEM 矩阵[14]。As-rigid-as-possible (ARAP) 是一种来自几何处理领域的网格变形方法，根据选定的变形点保持每个网格元素的质量[15]。ARAP 和其他几何处理方法适用于在接收到少数控制点后需要实时网格变形的应用[15]、[16]、[18]。然而，这些方法解决的基本问题是与 FEMWARP 不同的，因为变形后的边界表面是事先未知的。相比之下，FEMWARP 假设整个输入和输出边界几何是已知的，这对于前面提到的科学和工程应用非常重要。总体而言，本文讨论的网格变形类型是边界驱动的，而不是由控制点驱动的。我们提出了一种新颖的、仅基于 MPI 的 FEMWARP 并行化版本（ParFEMWARP）[1,2]，通过共享内存（SHM）实现节点内部通信，以及通过远程内存访问（RMA）实现节点间单边通信[19]。

其他网格变形方法旨在在应用大变形时生成高质量网格。这些方法基于非线性弹性建模变形[20]或执行网格细化以提高变形网格的质量[21]、[22]。这些方法在计算上更为昂贵，但可以轻松实现并行化。

本文是我们 2025 年 SIAM 国际网格研讨会上初步工作的扩展和邀请提交[23]。我们扩展了论文，提供了关于 ParFEMWARP 的更深入的技术和理论见解。具体来说，我们描述了分离物理和运动逻辑的抽象，以支持提供给 ParFEMWARP 的用户定义变形。此外，我们提供了组成 ParFEMWARP 每个阶段的算法的详细说明。我们还为方法的每个阶段提供了详细的并行运行时间和空间复杂度分析。这些补充提供了与 ParFEMWARP 在共享和分布式内存架构上变形大型四面体网格的能力相关的理论和实践细节。我们还扩展了文献综述，为我们的工作提供了额外的背景信息。

在本文中，我们描述了用于四面体网格的 ParFEMWARP。第 2 节描述了 FEMWARP 算法。第 3 节提供了有关 MPI-3 标准中混合并行编程的详细信息。第 4 节描述了我们如何利用 MPI-3 中可用的共享内存和单边通信特性来并行预计算邻居列表，这对于执行 FEMWARP 的各个阶段是必要的。第 5 节概述了我们用于解决多个右侧问题的预处理块共轭梯度方法。第 6 节描述了我们如何并行化 FEMWARP 的各个阶段。第 7 节提供了有关 ParFEMWARP 如何通过分离物理/运动逻辑与并行网格变形逻辑的抽象来支持用户定义变形的详细信息。我们在第 8 节提供了我们方法的并行运行时间和空间复杂度分析。第 9 节给出了我们的数值实验结果。第 10 节描述了结论和未来工作的方向。

**节选**
FEMWARP [1]、[7]、[8]、[9]、[24] 是一种基于有限元的网格变形算法。FEMWARP 遵循三步过程来计算变形后的网格。
FEMWARP 首先量化初始网格的表示。为此，它使用标准 Galerkin FEM 在域 Ω 上解决以下泊松边值问题 Δu=0 在 Ω上，u=u0 在 ?Ω上。可以根据 [7]、[8] 选择任何 u0。离散化产生全局刚度矩阵，其条目为 A(i,j)=∫Ω??i???jdΩ，其中 ?i 是分段线性函数。

**混合并行编程与 MPI-3**
在现代并行科学计算中，当在集群上并行化软件时，通常将 MPI 和 OpenMP 分别用于节点间和节点内通信。这种配对是混合并行编程的一个例子，它在节点和核心级别都实现了并行性。MPI-3 标准通过将 SHM 作为 MPI RMA 接口的一部分，提供了一种更统一的混合并行性替代方案，如图 1 所示。此外，Hoefler 等人展示了……

**预计算全局刚度矩阵的稀疏性**
矩阵 AI 以压缩稀疏行（CSR）格式存储，因为 AI 是一个大型稀疏矩阵，且是对称正定的（SPD）。向 CSR 矩阵插入新的非零元素从而修改其稀疏模式是非常昂贵的，并被强烈建议避免[27]。考虑到大型四面体网格的 FEM 全局刚度矩阵中可能有数十亿个非零值，因此需要预先计算 AI 的稀疏模式。

**预处理块共轭梯度**
预处理块共轭梯度（PCG）方法是一种用于解决大型、稀疏 SPD 线性系统的迭代方法[28]。然而，如果在 FEMWARP 中使用 PCG 解决多个右侧问题，则需要分别解决每个右侧的线性系统（即 AIx?I=?ABx?B, AIy?I=?ABy?B, 和 AIz?I=?ABz?B）。而预处理块共轭梯度（PBCG）方法[29]、[30] 可以用来解决多个右侧问题。

**我们的实现**
我们的实现包括第 4 节中描述的邻居预计算步骤，然后是第 2 节中显示的三步算法。注意我们没有使用网格划分器。如图 9 和算法 3 所示，基于等级局部四面体计算刚度矩阵贡献的过程是高度并行的。每个等级首先计算 AI 和 AB 的局部、部分版本（见公式 (2)–(3)）。请注意，不属于特定等级的边的条目将保持为零。

**用户定义的变形**
ParFEMWARP 支持来自各种来源的边界变形，例如来自实验数据、实时流式数据或指定变形的偏微分方程（PDE）。在算法层面，这是由于原始网格中的顶点连接性在整个变形过程中得以保留（即网格拓扑固定）。这带来了一个关键优势，即不需要进行昂贵的重新网格化计算。具体来说，与并行网格变形相关的计算开销较低。

**并行运行时间和空间复杂度分析**
在本节中，我们分析 ParFEMWARP，以提供我们并行化各个阶段的运行时间和空间复杂度的上限。设 t 表示四面体元素的数量，v 表示顶点的数量，d 表示任何给定顶点的最大度数，n 表示节点的数量，pglobal 表示均匀分布在 n 个节点上的进程数量，pnodal 表示节点内的进程数量，即 pnodal=pglobal/n，Cx 表示事件 x 发生所需的时间。假设 tt、tv、td 和 tn 分别表示对应的时间。

**数值实验**
我们通过运行一系列数值实验来评估 ParFEMWARP 的性能。ParFEMWARP 是使用 GNU C++、OpenMPI 4.0.7 和 Eigen 3.4.0 [32] 开发的。实验使用了 KU Community Cluster Bigjay 分区上的最多 256 个核心，这些核心配备有运行在 2.60 GHz 的 Intel Xeon Gold 6349 CPU 和 FDR Infiniband。节点每个包含 256 GB 或 1024 GB 的 RAM。我们使用两种网格进行实验：（1）手动制作的（图 11），包含大约 1 亿个四面体。

**结论和未来工作**
我们提出了 ParFEMWARP，这是一种用于变形四面体网格的 FEMWARP 算法的并行版本[7]、[8]。我们的方法利用了 MPI-3 中可用的共享内存和单边通信特性，因此实现可以在共享和分布式内存架构上移植。我们利用混合并行性来预计算全局刚度矩阵的稀疏性。具体来说，我们使用共享内存和单边通信来在计算邻居列表时获得良好的性能。

**作者贡献声明**
Abir Haque：写作——审稿与编辑、写作——原始草稿、可视化、验证、软件、方法论、调查、形式分析、数据管理、概念化。
Suzanne M. Shontz：写作——审稿与编辑、监督、资源获取、项目管理、方法论、资金获取、概念化。

**利益冲突声明**
作者声明以下可能被视为潜在利益冲突的财务利益/个人关系：Abir Haque 报告称国家科学基金会提供了财务支持。Suzanne M. Shontz 报告称国家科学基金会提供了支持。如果还有其他作者，他们声明没有已知的可能会影响本文所述工作的财务利益或个人关系。

**致谢**
这项工作是在堪萨斯大学（KU）研究中心运营的高性能计算（HPC）设施中进行的，该设施得到了国家科学基金会拨款 OAC-2117449 的支持。第一作者在隶属于堪萨斯大学电气工程与计算机科学系和信息科学研究所期间完成了大部分工作。第一作者的工作得到了 NSF 授款 OAC-1808553 的 REU 补助金支持。

**致谢**