《Journal of Computational Physics》：A predictor-corrector scheme for approximating signed distances using finite element methods
编辑推荐：
在本文中，研究人员介绍了一种有限元方法，旨在针对二维和三维空间中任意边界，稳健地计算近似符号距离函数。该方法采用了一种新颖的预测-校正方法，首先求解一个基于线性扩散的预测问题，随后求解一个与程函方程相关的基于非线性最小化的校正问题。预测步骤能够高效地生成合适
在本文中，研究人员介绍了一种有限元方法，旨在针对二维和三维空间中任意边界，稳健地计算近似符号距离函数。该方法采用了一种新颖的预测-校正方法，首先求解一个基于线性扩散的预测问题，随后求解一个与程函方程相关的基于非线性最小化的校正问题。预测步骤能够高效地生成合适的初始猜测值，显著促进了非线性校正步骤的收敛。研究人员使用具有任意陡峭或平坦区域的复杂初始水平集函数（在二维和三维空间中），以及瞬态测试案例（如涡旋测试），验证了该方法的稳健性。研究人员将结果与标准有限差分Hamilton-Jacobi格式进行了比较，证明了该方法在界面保持方面的改进。此外，研究人员还利用分布式内存并行（MPI）技术展示了该方法在大规模三维问题上的可扩展性，确立了其在重初始化各种水平集函数方面的准确性、高效性和广泛适用性。
**论文解读**
**一、研究背景与问题陈述**
符号距离函数（Signed Distance Function, SDF）的计算在计算力学、计算机图形学以及涉及界面捕获与追踪的数值模拟中扮演着至关重要的角色。其在形状优化、多相流模拟和图像处理等领域有着广泛的应用需求。符号距离函数通常被定义为程函方程（Eikonal Equation）的解，这是一个非线性偏微分方程，其数值求解面临重大的计算挑战。
现有的符号距离函数计算方法主要分为两大类：几何算法和基于偏微分方程（PDE）的方法。几何算法包括快速行进法（Fast Marching Method）和快速扫描法（Fast Sweeping Method），这些方法在结构化网格上高效且稳健，但难以并行化，且不易扩展到非结构化网格或高阶几何表示。基于PDE的方法，例如Sussman等人提出的时间依赖Hamilton-Jacobi（HJ）公式，通过在有限元、有限体积或有限差分框架下直接求解程函方程的变体，为复杂域提供了更大的灵活性，但需要仔细处理以保证稳定性和准确性。然而，当使用有限元方法（FEM）求解双曲型的HJ方程时，标准Galerkin格式不稳定，需要引入如流线迎风Petrov-Galerkin（SUPG）等稳定化技术，但这些技术会引入人工扩散，导致重初始化过程中出现界面位移和显著的质量损失。
为克服这些困难，近年来出现了基于最小化的重初始化框架。该类方法将重初始化问题转化为最小化梯度范数与1的偏差的最小二乘问题。尽管有效，但原始公式在处理梯度范数消失（平坦区域）的区域时会遇到奇异性，导致迭代不收敛。以往的研究通过修改扩散系数或势函数来正则化问题，但这些修改可能会改变被求解方程的性质，或无法完全解决任意初始条件下的收敛问题。因此，亟需一种既能保持原始一致的最小化问题，又能确保在任意初始条件下稳健收敛的新方法。
**二、研究方法与技术**
本研究所提出的预测-校正（Predictor-Corrector, P.C.）方法，其核心创新在于为解决非线性最小化问题的收敛困难提供了一个稳健的初始化步骤，而非修改程函方程本身。该方法主要包含两个关键步骤：
1. **预测步骤：** 构建一个线性扩散问题，其右端项和Neumann边界条件均基于符号函数（sign function）。该预测器问题的解析解特性保证了其能有效消除初始水平集函数中的平坦和陡峭区域，为后续的校正步骤提供一个具有非零梯度范数的良好初始猜测值。该设计源于对一维问题的解析推导。
2. **校正步骤：** 使用原始的、未修改的基于最小二乘最小化的非线性扩散问题（即Picard迭代格式）进行求解。该步骤通过迭代最小化梯度范数与1的偏差，逐步恢复符号距离性质。在校正步骤中，为处理梯度范数接近零的等距点，引入了一个极小值ε（约10?1?至10?12）作为数值稳定常数，以避免除零错误。
在离散化方面，研究采用了有限元方法。对于网格与界面不吻合的情况，研究采用了CutFEM方法。该方法通过Nitsche方法弱施加Dirichlet边界条件，并采用混合维度的子积分技术对切割单元和界面进行精确积分，从而避免了重新划分网格的需求。所有数值实现基于CutFEMx库（扩展自FEniCSx），并利用MPI实现分布式内存并行计算。
**三、研究结果**
研究通过一系列数值实验验证了该方法的性能，包括二维和三维的圆、阶跃函数、带噪声的界面、星形域以及三维圆环等形状。
* **经典算例：圆形水平集（Convergence study）**：通过对圆形界面函数的研究，确定了该方法的最优收敛阶。研究表明，经过约3次校正迭代，L2误差和程函误差即可降至较低水平（低于0.01和0.05），显示出极高的效率。在界面处的收敛阶为1.99，接近最优的2阶。在去除中心奇点（等距点）的影响后，所有误差度量的收敛阶均恢复至理论最优值。这表明，整体收敛阶的降低主要源于等距点处的不可微性，而非方法本身的缺陷。
* **Nitsche参数敏感性分析**：对Nitsche惩罚参数γ_D的分析表明，最终收敛时的程函误差与γ_D无关，证实了方法的稳健性。然而，较大的γ_D会降低收敛速度，但能更有效地保持质量守恒。
* **对比性能分析**：与标准的有限差分Hamilton-Jacobi（FD-HJ）求解器对比，P.C.方法虽然在总计算时间上更长，但其在界面保真度上表现优异，界面误差‖?‖_L2(Γ)约为1.1×10??，比最小二乘方法好近60倍，比HJ方法好4倍。这表明P.C.方法在界面保持方面具有显著优势，这对于水平集应用至关重要。
* **平坦函数与不规则函数测试**：面对梯度范数为零的大片平坦区域（如阶跃函数），无预测器的最小化方法会因右端项无穷大而发散。P.C.方法通过预测步骤成功生成了梯度范数非零的初始函数，从而使得校正步骤能够顺利收敛。对于带有随机扰动的反正切函数，P.C.方法显著优于直接使用最小化方法，能够平滑界面附近的显著不规则性，并随迭代次数增加不断改善误差。
* **复杂界面：星形水平集**：在具有高曲率和大量等距点的星形域测试中，P.C.方法经过少量迭代即可使程函误差显著下降。Nitsche方法确保了界面在重初始化过程中无可见位移，总体界面位移误差保持在O(10?3)量级。质量守恒通过预测步骤的微小外移和校正步骤的等量内移得以补偿实现。
* **三维模拟与并行性能**：在包含超过两百万未知量的三维圆环算例中，P.C.方法在6个处理器上仅用7分54秒即可得到高精度的符号距离函数近似。在严格的停滞准则（ε=10??）下，该方法的扩展性测试表明，其经验计算复杂度约为O(DoF1·1?)，接近理论上的线性复杂度，展示了其在大规模三维问题上的可行性和优越的并行可扩展性。
* **瞬态测试：单涡旋测试**：在Rider-Kothe涡旋测试中，P.C.方法在经历强烈的界面拉伸和拓扑变形后，仍能有效恢复初始界面形状。研究发现在强瞬态变形下，采用分段常数法向会导致界面附近产生振荡，通过引入基于L2投影的连续分段线性法向，可显著提高重初始化过程的稳定性。在最终时刻，非结构化网格的形状误差约为结构化网格的一半，显示出其各向同性特性在复杂变形中的优势。尽管质量误差和程函误差在最大拉伸时刻（t=1.5）达到峰值，但随着速度场反转，这些误差对称地恢复到初始水平，表明该方法在可逆变形下不会积累数值噪声。
**四、结论与讨论**
本研究提出了一种基于有限元方法的预测-校正格式，用于从任意形状的水平集函数稳健地计算符号距离函数。该方法通过一个线性扩散预测器和一个非线性最小化校正器，有效解决了传统最小化方法在梯度消失区域的收敛问题，同时避免了求解双曲型HJ方程所需的复杂稳定化技术。研究结论可总结如下：
1. 该方法对包括光滑、不连续和不规则在内的多种初始水平集函数均表现出稳健性和收敛性，并能有效处理三维复杂几何形状。
2. 该方法易于并行化，在处理数百万未知量的大规模网格时，仅需少量迭代即可完成水平集函数的重初始化。
3. 通过涡旋测试验证了该方法在瞬态问题中的稳健性，确认了其在可逆变形下能有效保持界面几何形状并确保质量守恒。
该研究的核心贡献在于提出了一种不修改原始最小化问题，而是通过精心设计的预测步骤提供稳健初始猜测的策略，从而在保持问题一致性的同时解决了收敛难题。其基于CutFEM的框架天然支持非结构化网格和并行计算，为大规模、复杂几何形状下的水平集重初始化提供了一种高效且鲁棒的替代方案，尤其在需要高水平界面保真度的应用场景中具有重要价值。
