《ARCHIVES OF COMPUTATIONAL METHODS IN ENGINEERING》：Comparative Study of Different Quadrature Methods for Cut Elements
切割单元求积对于所有不采用边界拟合网格的有限元方法（Finite Element Method，FEM）至关重要。它应当高效、精确且鲁棒。已有多种平衡这些需求的方法被发表，其中一些以开源实现形式提供。本研究回顾了这些开源代码及其所使用的方法。此外，研究人员还针对二维和三维几何体开发了基准测试示例。所有模型均支持隐式和显式边界描述。不同的示例测试了代码的效率、精度、通用性和鲁棒性。特别关注了控制期望求积阶（quadrature order）的输入参数对实际积分误差的影响。对讨论的代码进行了详细比较。系统性的基准测试使得能够进行结论性比较，并为未来的代码开发提供了有价值的工具。所有测试均发布在随附的开源仓库中。
**1 引言（Introduction）**
有限元方法（Finite Element Method，FEM）在研究和工程中广泛用于各种模拟，其预处理阶段需要对几何体进行鲁棒的网格剖分，但对任意复杂几何体而言，网格剖分可能成为繁琐且耗时的任务，因为常需手动干预（如几何修复、去特征和重复网格剖分）。浸入式方法（immersed methods）旨在规避这一操作，其基本思想是将（复杂）几何体嵌入一个简单的背景网格中，避免上述网格操作。控制方程在较简单的均匀笛卡尔背景网格上离散，网格被域边界分割为活跃区域和非活跃区域，部分单元可能被域边界穿过——这类单元称为切割单元（cut elements），因此需要特殊处理。众多方法实现了这一基本思想的变体，如CutFEM、嵌入域方法、扩展有限元法（extended FEM，XFEM）、虚拟域方法、有限胞元法（Finite Cell Method，FCM）、非适配有限元等。这些方法共享共同的发展历史：连续扩展有限元法起源于固体力学中的裂纹问题（Mo?s等人），之后推广至时间相关问题、两相流问题；内在XFEM公式（Fries）实现了速度场中的折线表示；扩展间断伽辽金（extended Discontinuous Galerkin，XDG）方法（又称非适配或切割胞元DG方法）由Bastian和Engwer提出，后扩展至两相流；CutFEM框架（Burman等人）应用于涉及复杂和演变形体的广泛问题；有限胞元法（Parvizian等人）结合了虚拟域思想与高阶近似；等几何分析（Isogeometric Analysis，IGA）（Hughes等人）统一了基于CAD的几何表示和有限元分析。浸入式方法的概念也出现在许多通过修剪操作构建的CAD模型中。尽管所列方法在名称或公式细节上可能不同，但本研究广泛采用“切割单元”这一术语以突出它们在此类单元上的共同问题和任务。浸入式公式简化了模型准备，但引入了与标准边界拟合FEM不同的数值任务。综述[19]列出了三个核心挑战：切割单元上积分数值评估、浸入式或非适配边界上（本质）边界条件的施加、与小切割单元相关的格式稳定性。本研究处理第一个挑战——切割单元上的求积。背景网格的张量积结构被任意切割边界破坏，因此需要在被切割单元上应用定制的求积规则。大量研究致力于此任务，可能的策略包括四叉树/八叉树细分、矩拟合、单元重新参数化、多项式积分或积分降维。众多方法使选择合适工具变得复杂，所幸部分开发的方法已作为开源代码发布，可直接访问并快速应用。本文回顾了可用的代码，重点关注它们用于切割单元求积的方法，共找到13个开源仓库，其中10个被详细讨论和测试。研究者实现了一个Matlab基准环境，可从共同基础执行这些代码，该接口已开源。此前的一篇论文[24]讨论了接口结构，展示了基于持续集成（Continuous Integration，CI）的工作流用于测试和基准测试，并从软件工程角度开发了基准测试基础设施。本文则侧重于系统性的基准选择，以全面比较不同代码及其应用方法，重点测试其精度、鲁棒性、通用性和效率。统一符号便于不同方法之间的比较。在浸入式方法中，边界通常由隐式或参数化表示定义，因此研究者设计的基准可被两种表示精确定义，从而允许比较使用其中任一表示的工具。水平集函数（level set functions）被视为隐式表示，B样条（B-spline）或非均匀有理B样条（Non-Uniform Rational B-Spline，NURBS）作为参数化描述。基准测试考虑均匀笛卡尔背景网格，涵盖2D和3D示例。这些基准不仅使现有方法能够客观比较，也为研究软件中“测试案例设计”这一广泛认可的挑战提供了标准化测试集，用于评估未来实现。
**2 问题描述与符号（Problem Description and Notation）**
主要关注切割单元的求积，这些单元定义任意域\(K \subset \mathbb{R}^{{d^e}}\)（在\(\mathbb{R}^d\)中）且具有光滑参数化\(F: K \rightarrow \mathbb{R}^d\)，其中\(d^e \leq d\)，\(d^e\)是单元维数，\(d\)是空间维数。本节回顾基于文献[28-33]的逼近理论，首先介绍数值求积公式的基本符号和概念，然后总结域变换如何影响积分精度（尤其是需要几何近似时），为评估切割单元求积方法的精度提供理论基础，最后讨论所需求积点数的确定。
**2.1 基本求积问题（The Basic Quadrature Problem）**
“求积”（quadrature）指数值计算积分。最简单形式是计算参考域\(\hat{K}\)上连续函数\(\hat{f}\)的积分\(I(\hat{f}) = \int_{\hat{K}} \hat{f}(\hat{\boldsymbol{x}}) d\hat{\boldsymbol{x}}\)。求积方法通过一组求积点\(\{\hat{\boldsymbol{x}}_1, ..., \hat{\boldsymbol{x}}_{n_{QP}}\}\)和对应权重\(\{\hat{w}_1, ..., \hat{w}_{n_{QP}}\}\)来逼近\(I\)。求积阶\(k_q\)是多项式空间\(\mathbb{P}_k\)中使公式精确成立的最大次数\(k\)。误差可表示为\(|\hat{E}_{k_q}| \leq c h_{\hat{K}}^{k_q+1} |\hat{f}|_{k_q+1,\infty,\hat{K}}\)，其中\(h_{\hat{K}}\)是\(\hat{K}\)的直径，\(|\cdot|_{k_q+1,\infty,\hat{K}}\)是Sobolev空间\(W^{k_q+1,\infty}(\hat{K})\)中的半范数。
**2.2 域变换（Domain Transformations）**
设\(K\)由映射\(\mathbf{T}_K: \hat{K} \rightarrow K\)定义，积分换元后得到\(I(f) = \int_K f(\boldsymbol{x}) d\boldsymbol{x} = \int_{\hat{K}} f(\mathbf{T}_K(\hat{\boldsymbol{x}})) \det(\nabla\mathbf{T}_K) d\hat{\boldsymbol{x}}\)。对于切割单元，通常用近似\(\tilde{K}\)（由\(\tilde{\mathbf{T}}_K: \hat{K} \rightarrow \tilde{K}\)给出）代替真实域\(K\)。总积分误差分解为几何近似误差\(E_T\)和数值积分误差\(E_Q\)。若\(\tilde{K}\)由\(\mathbf{T}_K\)的多项式插值（次数\(q\)）指定，则\(|E_T| \leq c h_{\hat{K}}^{q+1} |\mathbf{T}_K|_{q+1,\infty,\hat{K}}\)；而\(|E_Q| \leq c h_{\hat{K}}^{k_q+1} |\tilde{\gamma}|_{k_q+1,\infty,\hat{K}}\)。总体误差为\(\left|I(f) - \tilde{I}_n(f)\right| \leq c h_{\hat{K}}^{\kappa}\)，其中\(\kappa = \min\{k_q+1, q+1\}\)。对于非多项式边界，\(E_T\)构成几何误差下限，而精确表示几何的方法可完全消除\(E_T\)。
**2.3 求积点数（Number of Quadrature Points）**
**2.3.1 求积公式的求积阶**：求积阶\(k_q\)与求积点数\(n_{QP}\)有关，例如Gauss公式达到\(k_q=2n_{QP}-1\)。
**2.3.2 域变换的影响**：当存在域变换时，被积函数\(f \circ \tilde{\mathbf{T}}_K\)的代数次数增加。对于一般映射、带一条曲边的四边形映射、轴对齐带一条曲边的四边形映射、退化映射以及仿射变换，给出了被积函数与雅可比行列式乘积的张量积多项式次数。这些次数直接影响所需Gauss求积点数。关键结论是：几何映射的多项式次数\(p\)不足以准确估计主导次数\(d_{\max}\)，而知道\(d_{\max}\)对设置正确求积点数至关重要。对于张量积Gauss规则，要求\(2n_{QP}-1 \geq d_{\max}\)，即\(n_{QP} = \lceil (d_{\max}+1)/2 \rceil\)。
**3 开源代码综述（Review of Open-Source Codes）**
本节介绍所有考虑的开源代码。首先在**3.1 概述**中给出13个开源代码的概览，按切割单元求积方法分组，包括重新参数化、矩拟合、降维等类别。大多数代码处理隐式边界表示，少数支持参数化描述。所有代码可处理2D和3D几何，但QUGaR和QuaHOG除外。应用领域和编程语言多样。
**3.2 三个2D测试几何**：引入三个不同曲线次数（线性、二次、五次）的2D几何，用于后续单代码研究。它们模拟单个单元的切割情况，面积可解析计算。研究考察三个问题：增加求积点数\(n_{QP,set}\)对面积计算精度的影响、\(h\)-细化对面积计算精度的影响、增加\(n_{QP,set}\)对单项式积分精度的影响。对于高阶求积方案（Algoim、BoSSS、ngsxfem、QuaHOG、QUGaR），最优\(n_{QP,set}\)取决于域变换，轴对齐四边形映射下所需点数按公式\(n_{QP,set} = \lceil (qp+q+p+1)/2 \rceil\)估算。
**3.3 四叉树/八叉树（Quadtree/Octree）**：该方法通过递归细分切割单元并在子单元中放置高斯点，配合内部/外部检查。
**3.3.1 FCMLab**：应用四叉树/八叉树方法，使用隐式边界，两种边界积分方法（低精度使用四叉树，高精度使用线性边界近似）。关键参数：高斯点数、细分层数（本研究用3层）、线性边界近似种子点数。研究显示，对于三种几何，面积计算误差随\(n_{QP,set}\)增大对情况1和2降低，但情况3无变化；\(h\)-细化收敛约二阶；单项式积分精度随\(n_{QP,set}\)提高，但远未达机器精度。
**3.3.2 Nutils**：应用四叉树/八叉树方法，并在最低细分层上附加三角剖分，使用单纯形求积规则，仅支持隐式边界。关键参数：求积阶、细分层数（本研究用3层）。面积计算误差不随\(n_{QP,set}\)变化，但随边界曲线次数降低而降低（因三角剖分固定）；\(h\)-细化呈现稳定二阶收敛；单项式积分精度随\(n_{QP,set}\)提高，但受几何误差限制。
**3.4 重新参数化（Reparametrization）**：此类方法通过映射标准求积规则到参数化良好的切割单元。
**3.4.1 Gridap**：对隐式边界进行线性近似，在单元角点处评估水平集函数并线性插值，然后对线性边界的切割单元进行三角剖分，使用映射高斯积分（2D三角形，3D四面体）。关键参数：求积阶。无法用单个单元解析几何，仅进行\(h\)-细化研究，面积计算呈现二阶收敛（线性边界可达机器精度）。
**3.4.2 ngsxfem**：与Gridap类似，但通过高阶网格变换改进边界近似。优选三角形/四面体网格，也支持四边形网格但限于单水平集函数。关键参数：求积阶、重新参数化次数（本研究设为曲线次数）。线性边界可精确解析，对二次和五次边界进行\(h\)-细化，可达高阶收敛（但重新参数化次数需匹配曲线复杂度）；单项式积分在情形1中可达机器精度。
**3.4.3 QUGaR**：独立C++库，主要用于参数化界面，但隐式边界也可用。通过模板将切割单元分为四边形子区域并建立直纹曲面，必要时分裂单元。关键参数：高斯点数、重新参数化次数。三种几何的面积和单项式积分均可在预期\(n_{QP,set}\)下达到机器精度（五边形单元需分裂导致双倍点数）。
**3.5 矩拟合（Moment Fitting）**：通过求解线性方程组确定求积权重，节点位置固定。
**3.5.1 BoSSS**：基于扩展间断伽辽金方法，采用分层矩拟合（hierarchical moment fitting，HMF）。利用散度定理将体积积分转化为边界积分，再通过矩拟合求解。引入安全系数\(\varpi\)控制节点数。关键参数：求积阶、安全系数。面积计算可在\(n_{QP,set}\geq2\)时达机器精度（但实际节点数较多）；单项式积分同样收敛，但对于非多项式被积函数（如情形3）无法达机器精度，误差低于\(10^{-5}\)。
**3.5.2 mlhp**：FCMLab的继任者，支持hp-有限元，采用矩拟合方法克服四叉树/八叉树点数过多的问题。节点选为Gauss-Legendre点，基函数为Lagrange多项式，使矩阵成为单位矩阵。关键参数：最终矩拟合高斯点数、四叉树参考积分高斯点数、细分层数、种子点数。面积计算误差与FCMLab几乎一致（因使用四叉树参考），但\(n_{QP}\)显著减少；\(h\)-细化收敛约二阶；单项式积分精度随\(n_{QP,set}\)提高，但远未达机器精度。矩拟合的准确性依赖于参考积分的精度。
**3.5.3 QuESo**：针对3D嵌入几何的预处理器，接收STL文件。使用矩拟合方法，基函数为3D Legendre多项式，通过点消除算法获得最优求积点。利用散度定理计算精确积分，边界由STL三角形网格表示。关键参数：求积阶、矩拟合残差容差。仅支持3D，因此未进行2D几何研究。
**3.6 降维（Dimension-Reduction）**
**3.6.1 Algoim**：基于高度函数降维，将积分转化为嵌套积分，递归应用一维Gauss-Legendre规则。隐式边界定义，通过水平集函数梯度方向确定高度方向。非单调时需细分。关键参数：求积点数、重新参数化次数。三种几何的面积和单项式积分均可在预期\(n_{QP,set}\)下达到机器精度（五边形单元需细分导致双倍点数）。
**3.6.2 QuaHOG**：基于Green定理将二维积分降为一维线积分，边界由有理Bernstein-Bézier曲线组成。通过中间求积和反导数求积两步得到最终求积点。关键参数：中间求积点数\(Q_i\)、反导数求积点数\(P\)。默认设\(Q_i=P\)。附加变体QuaHOGPE使用有理求积规则，对任意多项式积分精确。面积和单项式积分在预期\(n_{QP,set}\)下达到机器精度，QuaHOGPE甚至对NURBS曲线也达到机器精度。
**3.7 其他代码**：CutFEM-Library、deal.II和TPMC未包含在详细比较中。CutFEM-Library和deal.II直接使用Algoim；TPMC提供拓扑保持移动立方体算法，但缺少充分文档和示例。
**3.8 计算效率**：由于编程语言、接口开销等因素，无法定量比较计算时间，仅定性讨论预处理成本。四叉树/八叉树方法只需内部/外部检查，成本低但点数多；三角剖分和网格变换增加预处理时间；矩拟合需求解线性系统并计算参考积分，成本最高；降维方法（Algoim）需确定高度函数；QuaHOG直接基于曲线积分，成本较低。
**4 比较性数值示例（Comparative Numerical Examples）**
本节通过多个数值示例评估和比较不同开源源代码及其求积方法，包括纯数值积分（无特定物理应用），分析几何误差和积分误差，涵盖2D和3D几何。通用说明：FCMLab与mlhp结果几乎一致；Gridap无法处理单元素水平集；ngsxfem支持四边形网格和网格变换（仅单水平集），3D仅四面体；Nutils无法区分切割与未切割单元的求积点；QuaHOG和QuaHOGPE不参与\(h\)-细化；QuaHOG中设\(Q_i=P=n_{QP,set}\)。
**4.1 二维几何（Two-Dimensional Geometries）**
**4.1.1 求积点布局**：展示情形3的求积点分布。Algoim和QUGaR在五边形单元中需要细分，FCMLab和Nutils显示四叉树高密度点，BoSSS、Gridap、mlhp、QuaHOG和QuaHOGPE的点可能位于活跃域外。
**4.1.2 多样化切割场景**：测试26种不同几何（不同曲线次数、有理曲线、尖点、内部节点、小特征等），进行面积计算（\(n_{ele}=1\)和\(n_{ele}=8\)）。代码分为两组：存在几何近似误差的（FCMLab、Gridap、mlhp、ngsxfem、Nutils）和精确表示边界的（Algoim、BoSSS、QUGaR、QuaHOG、QuaHOGPE）。多数代码在精细网格下相对误差低于\(10^{-5}\)，QuaHOGPE对所有26种几何达到完美精度。
**4.1.3 平移交会抛物线的鲁棒性测试**：将两个交会抛物线平移1000步，面积计算误差显示所有代码鲁棒，QUGaR达到机器精度且无振荡。统计均值误差和标准差，Gridap和ngsxfem误差较大，mlhp较FCMLab大幅减少求积点数，QUGaR在精度和点数上表现突出。
**4.1.4 单项式作为被积函数**：针对情形2（二次曲线）和六次单项式，研究\(n_{QP,set}\)和\(h\)-细化。区分几何误差主导的代码（FCMLab、mlhp、ngsxfem、Nutils）和准确处理几何的代码（Algoim、BoSSS、QUGaR、QuaHOG、QuaHOGPE）。后者在单元素下随\(n_{QP,set}\)提高精度，前者因几何误差限制无法受益。\(h\)-细化下，Gridap和Nutils二阶收敛，ngsxfem三阶（因网格变换），Algoim、BoSSS、QUGaR可达高阶收敛（\(2n_{QP,set}\)阶），FCMLab和mlhp二阶且受四叉树影响。
**4.2 三维几何（Three-Dimensional Geometries）**：支持3D的代码有Algoim、BoSSS、FCMLab、Gridap、mlhp、ngsxfem、Nutils和QuESo。
**4.2.1 平移椭球的鲁棒性测试**：将椭球平移1000步，体积计算误差显示所有代码鲁棒。Nutils在默认设置下893/1000步失败，需调整容差。QuESo不受平移影响（精度由STL弦公差决定）。统计均值误差表明Gridap仍存在几何分辨问题，mlhp较FCMLab平均减少约7.5倍求积点数。
**4.2.2 圆环的体积计算**：对NURBS描述的二次圆环进行\(h\)-细化。Gridap和Nutils二阶收敛，Algoim趋势接近四阶，BoSSS最后一步超出预期，ngsxfem三阶但波动，FCMLab和mlhp精度相同但点数不同，Nutils在精细网格上失败，QuESo不受网格细化影响。
**4.2.3 椭圆圆柱的单项式被积函数**：研究零到五次单项式，\(n_{QP,set}\)影响在3D中弱于2D，因几何误差占主导。仅展示\(q=3\)和\(q=5\)的结果，除Algoim外其他代码均因几何近似误差而饱和。\(h\)-细化下（\(q=5\)），Gridap二阶，BoSSS四阶，Algoim预渐近后趋于最低误差，ngsxfem三阶波动，FCMLab和mlhp一致。
**5 结论与展望（Conclusion and Outlook）**
开源代码提高了可重复性并加速研究。本文识别、回顾并比较了有关切割单元的可免费获取代码，涵盖多种切割单元数值求积方法，通过2D和3D域积分基准测试评估其通用性、鲁棒性、精度和效率。所有代码在测试场景中均表现鲁棒且可直接使用。特别研究了所需求积阶（即求积点数\(n_{QP,set}\)）对2D和3D示例的影响。提供了一个Matlab开源接口[23]，可从统一测试环境运行所有代码，并可检索大部分代码的最终求积点，便于在其他求解器中重用。回顾揭示出开源程序可分为两类：对边界进行三角剖分（引入几何误差，二阶收敛）和保留高阶边界或高阶近似（可达高阶收敛）。目前仍存在明显空白：缺乏结合参数化边界描述与体积的开源代码（尤其适用于B-Rep的CAD工作流）。未来可进一步研究改进参数（如重新参数化次数、细分层数、矩拟合容差），并比较浸入式方法中边界积分的不同处理。欢迎读者利用测试环境[23]对自己的实现进行基准测试。
