基于非均匀有理B样条（NURBS）几何的自适应数值积分策略
题目描述
考虑在计算几何与等几何分析中，物理量（如质量、刚度矩阵元素）常需在由非均匀有理B样条（NURBS）描述的曲线上或曲面域上进行积分。NURBS参数映射通常高度非线性，且被积函数在参数域中可能出现剧烈变化（如峰值、边界层），直接使用均匀网格的高斯求积效率低下。本题要求设计一种自适应数值积分策略，用于精确计算NURBS参数域上的积分，并能根据被积函数在物理域中的行为自动调整积分节点的分布。
具体问题：
设二维NURBS曲面由参数映射 \(\mathbf{x}(u,v) : [0,1]^2 \to \mathbb{R}^2\) 定义，需计算积分
\[I = \iint_{\Omega} f(\mathbf{x}) \, d\mathbf{x} = \int_0^1 \int_0^1 f(\mathbf{x}(u,v)) \, |J(u,v)| \, du\,dv, \]
其中 \(|J(u,v)|\) 是映射的雅可比行列式，\(f\) 是定义在物理域上的函数（可能包含峰值或振荡）。目标是设计自适应算法，在参数域 \([0,1]^2\) 上自动细分网格，并在子区域上应用高斯求积，使积分误差满足给定容差。
解题过程循序渐进讲解
1. 理解NURBS积分挑战
- 参数映射的非均匀性：NURBS的 \(u,v\) 参数到物理坐标的映射非线性较强，均匀参数步长可能对应物理域中扭曲或疏密不均的网格。
- 雅可比行列式的影响：\(|J(u,v)|\) 可能变化剧烈（例如在曲面高度弯曲处），导致被积函数 \(g(u,v) = f(\mathbf{x}(u,v)) |J(u,v)|\) 在参数域中表现出奇异性或剧烈振荡。
- 直接均匀求积低效：若在整个参数域上用固定高斯节点，可能需要极多节点才能捕获细节，计算量过大。
2. 自适应策略框架
核心思想：在参数域上进行递归区域细分，基于误差估计决定是否细分。步骤包括：
- 初始网格：将参数域 \([0,1]^2\) 划分为若干粗网格（如基于NURBS节点向量生成初始细分）。
- 子区域积分估算：在每个子矩形区域上应用二维高斯求积（如 \(m \times m\) 点高斯-勒让德求积）计算积分近似值。
- 误差估计：通过比较不同精度的求积结果（例如低阶与高阶高斯求积）得到局部误差估计。
- 细分判据：若子区域的误差超过允许容差（考虑区域面积加权），则将该区域划分为四个相等子矩形。
- 递归执行：对新区域重复步骤2-4，直至所有区域满足误差要求或达到最大细分深度。
- 汇总结果：将所有子区域的积分值求和，作为总积分近似。
3. 关键技术细节
（1）NURBS参数映射与雅可比计算
- 给定NURBS的控制点、权重、节点向量、基函数阶数，可计算任意 \((u,v)\) 处的坐标 \(\mathbf{x}(u,v)\) 及雅可比矩阵 \(J(u,v) = \partial \mathbf{x} / \partial (u,v)\)。
- 雅可比行列式 \(|J(u,v)|\) 需精确计算，因为它决定了积分测度变换。可通过B样条基函数的导数公式得到。
（2）二维高斯求积在参数子区域上的应用
- 对于参数域子矩形 \([a,b] \times [c,d]\)，通过变量变换将其映射到标准区间 \([-1,1]^2\)：
\[ u = \frac{b-a}{2} \xi + \frac{a+b}{2}, \quad v = \frac{d-c}{2} \eta + \frac{c+d}{2}. \]
- 积分近似为：
\[ I_{\text{sub}} \approx \frac{(b-a)(d-c)}{4} \sum_{i=1}^m \sum_{j=1}^m w_i w_j \, g\left(u(\xi_i), v(\eta_j)\right), \]
其中 \((\xi_i, w_i)\) 是 \(m\) 点高斯-勒让德求积的节点与权重，\(g(u,v)=f(\mathbf{x}(u,v))|J(u,v)|\)。
（3）局部误差估计的实用方法
- 常用嵌入误差估计：在同一子区域上分别用较低阶 \(m_1\) 和高阶 \(m_2\)（\(m_2 > m_1\)）高斯求积计算积分值 \(I_{\text{low}}\) 和 \(I_{\text{high}}\)，以两者差值作为误差估计：
\[ E_{\text{sub}} = |I_{\text{high}} - I_{\text{low}}|. \]
- 为平衡全局误差，将容差 \(\varepsilon\) 按子区域面积比例分配。设总容差为 \(\varepsilon_{\text{global}}\)，当前区域面积（参数域中）为 \(A_{\text{sub}}\)，总参数域面积为 1，则允许的局部误差为 \(\varepsilon_{\text{sub}} = \varepsilon_{\text{global}} \cdot (A_{\text{sub}} / 1)\)。
- 若 \(E_{\text{sub}} > \varepsilon_{\text{sub}}\)，则标记该区域需细分。
（4）自适应细分策略
- 细分时沿 \(u\) 和 \(v\) 方向各平分，产生四个子区域。
- 为避免无限细分，设置最大细分深度（如最多将参数区间细分 10 层），或最小区域边长限制。
- 细分后对新区域重新计算积分和误差，直至所有区域满足精度要求。
（5）处理被积函数的奇异性或峰值
- 在物理域中，\(f(\mathbf{x})\) 的剧烈变化会通过映射影响 \(g(u,v)\)。自适应策略可自动在参数域中对应位置加密节点。
- 若已知 \(f\) 在物理域中某些位置有奇异性，可在初始化时预先在这些位置对应的参数区域进行细分，加速收敛。
4. 算法伪代码
输入：NURBS映射 x(u,v)，函数 f，全局容差 eps，最大细分深度 max_depth
输出：积分近似值 I_total
初始化：区域列表 regions = { [0,1]×[0,1] }，I_total = 0
while regions 非空:
取出一个区域 R
在 R 上计算低阶积分 I_low 和高阶积分 I_high
计算局部误差 E = |I_high - I_low|
计算允许误差 eps_local = eps * (R的面积)
if E <= eps_local 或 R 深度 == max_depth:
I_total += I_high # 接受高阶结果
else:
将 R 细分为4个子区域，加入 regions
5. 注意事项与优化
- 缓存机制：在细分过程中，可缓存 NURBS 基函数及其导数在参数节点上的值，避免重复计算。
- 并行化：各子区域的积分计算相互独立，可并行处理以加速。
- 误差控制可靠性：嵌入式误差估计在高振荡或奇异函数中可能不可靠，可辅以相对误差控制或参考区间细分前后的变化。
6. 示例说明
假设有一NURBS曲面描述一个带尖锐曲率的几何，需积分物理域中的电场能量密度 \(f(\mathbf{x}) = \exp(-100\|\mathbf{x} - \mathbf{x}_0\|^2)\)（一个高峰值函数）。
- 参数域中，峰值位置对应特定 \((u_0,v_0)\)。自适应算法会在该处自动细分，在其他平缓区域用较粗网格，从而在保证精度的同时减少求积节点总数。
总结
本题通过结合 NURBS 参数映射的几何特性与自适应递归细分，实现了对物理域上复杂被积函数的高效数值积分。关键在于利用参数域上的误差估计驱动局部网格加密，使计算资源集中在被积函数变化剧烈的区域，从而平衡精度与计算成本。该方法可推广至三维NURBS体积分，是等几何分析中常用的积分策略。