Warp 多项式工具指南：`warp.fem.polynomial` 的一维求积规则与 Lagrange 缩放因子

Warp 多项式工具指南：warp.fem.polynomial的一维求积规则与 Lagrange 缩放因子

【免费下载链接】warpA Python framework for GPU-accelerated simulation, robotics, and machine learning.项目地址: https://gitcode.com/GitHub_Trending/warp/warp

本篇技术指南以 warp.fem.polynomial API 参考 为骨架，围绕其公开的两个核心接口quadrature_1d与lagrange_scales展开。它服务于 Warp 有限元（FEM）框架中最基础也最关键的两项任务：在参考单元上生成数值积分点与权重（Gauss–Legendre、Lobatto–Gauss–Legendre 与等距 Newton–Cotes 家族），以及为 Lagrange 插值形函数预计算缩放系数。读完本文，你将掌握warp.fem.polynomial的完整 API 语义、四种多项式族的具体数值规则、底层源码实现原理，以及它们如何被张量积化为一维/三维求积、并被形函数与RegularQuadrature等求积类实际消费。

模块定位：API 文档背后的实现链

docs/api_reference/warp_fem_polynomial.rst是一份指向warp.fem.polynomial模块的自动 API 引用页，它公开导出两个函数：

• lagrange_scales：Lagrange 多项式的缩放因子；
• quadrature_1d：一维求积点与权重。

模块本身是薄封装层：warp/fem/polynomial.py 仅做重新导出，真正的实现位于 warp/_src/fem/polynomial.py：

from warp._src.fem.polynomial import quadrature_1d as quadrature_1d from warp._src.fem.polynomial import lagrange_scales as lagrange_scales

从源码结构看，整个warp.fem的求积（quadrature）与形函数（shape function）体系都以这两个函数为基石：Polynomial枚举同时被 quadrature 框架、参考单元实现 以及 square/cube 形函数 导入。这一点在测试代码中体现得最直观——warp/tests/fem/test_fem_quadrature.py中直接以fem.Polynomial.GAUSS_LEGENDRE的形式访问该枚举。

Polynomial枚举：四种一维多项式族

warp.fem.polynomial的核心概念是"多项式族"（Polynomial family），它决定了一维区间上的插值节点/求积节点布局。源码 定义了四种：

辅助函数is_closed(family)判断某族是否包含区间端点，实现为：

def is_closed(family: Polynomial): return family == Polynomial.LOBATTO_GAUSS_LEGENDRE or family == Polynomial.EQUISPACED_CLOSED

这一判定直接决定形函数的节点分类：闭型族把节点放在单元顶点/边/面上（顶点处VERTEX_NODE_COUNT = 1），而开型族（如 Gauss–Legendre）所有节点都落在单元内部（见 square_shape_function.py 与 cube_shape_function.py 的节点计数逻辑）。

quadrature_1d：一维求积规则 API

quadrature_1d(point_count, family)返回(coords, weights)二元组：

• coords：NumPy 数组，求积点坐标，统一归一化到[0, 1]区间；
• weights：NumPy 数组，对应权重，按区间长度归一化（∫₀¹ 1 dx = 1）。

import warp.fem as fem import numpy as np # 3 点 Gauss–Legendre 规则 coords, weights = fem.polynomial.quadrature_1d(3, fem.Polynomial.GAUSS_LEGENDRE) print(coords) # [0.5, 0.11270167, 0.88729833] print(weights) # [0.44444444, 0.27777778, 0.27777778] print(np.sum(weights)) # 1.0，权重和恒等于区间长度

注意第一个参数是求积点个数（point_count），而不是多项式阶数。阶数与点数的换算由参考单元层负责（见下文"从阶数到点数"一节）。

Gauss–Legendre 族

_gauss_legendre_quadrature_1d硬编码了 n=1 到 n=5 的经典高斯节点，内部先在[-1, 1]上构造，再整体移位缩放至[0, 1]：

n 个 Gauss–Legendre 点对不超过 2n−1 阶多项式精确成立（n=3 时权重和为 1，可精确积分到 5 阶），是精度/点数比最高的规则。它的代价是不含端点——用于装配时需要在单元边界做额外的通量处理（例如 DG 方法）。

Lobatto–Gauss–Legendre 族

_lobatto_gauss_legendre_quadrature_1d支持 n=2 到 n=5，强制包含区间两端点：

LGL 的端点节点使其成为连续有限元（如 C⁰ Lagrange）的天然选择：节点与单元顶点重合，方便全局自由度编号与装配。其精度为 2n−3 阶（n=3 时对应 Simpson 法则的 3 阶精度），略低于同点数的 GL，但换来端点插值能力。

等距族（Newton–Cotes）

等距族使用均匀分布的节点。quadrature_1d对EQUISPACED_CLOSED/EQUISPACED_OPEN分别分发到_closed_newton_cotes_quadrature_1d与_open_newton_cotes_quadrature_1d，权重取自经典的 Newton–Cotes 公式（源码注释引用了 MathWorld 与 OEIS A093735/A093736）。

闭型代表规则（归一化后）：

• n=2：[0.5, 0.5]（梯形法则，1 阶精度）；
• n=3：[1/6, 2/3, 1/6]（Simpson 法则，3 阶精度）；
• n=4：[1/8, 3/8, 3/8, 1/8]（Simpson 3/8 法则）；
• n=5：[14/180, 64/180, 24/180, 64/180, 14/180]（Boole 法则）；
• 更高阶 n 直至 8 均有硬编码权重（如 n=8 时中心权重含343/640）。

开型族（节点避开端点）则出现负权重：n=3 时[2, −1, 2]/3、n=5 时[11, −14, 26, −14, 11]/20。负权重意味着积分结果对舍入误差更敏感，高次等距 Newton–Cotes 也因 Runge 现象而不适合高阶使用——这正是源码将高阶精度任务交给 Gauss 族的原因。需要指出的是，等距族在数值上不如高斯族稳定，实践中主要服务于等距网格与简单测试场景。

支持的点数范围

四种族的实现覆盖点数有限，超出即抛NotImplementedError：

• GL：1–5 点；
• LGL：2–5 点；
• 闭型 Newton–Cotes：2–8 点；
• 开型 Newton–Cotes：1–7 点。

从阶数到点数：_point_count_from_order的换算规则

quadrature_1d直接吃"点数"，而有限元语境通常从"多项式阶数"出发。参考单元层在 element.py 中提供了_point_count_from_order(order, family)完成换算：

例如order=2、family=GAUSS_LEGENDRE得到 2 个点（两点 Gauss 规则可精确积分 3 阶多项式，满足二次被积函数需要）；而order=2的 LGL 需要2//2+2 = 3个点（Simpson 法则）。这一映射保证了"积分精度不低于被积多项式的阶数"。

lagrange_scales：Lagrange 形函数的缩放因子

lagrange_scales(coords)接收一组节点坐标，返回对应的 Lagrange 基函数缩放系数。对第 i 个节点，其值定义为：

scale[i] = 1 / ∏_{j ≠ i} (coords[i] − coords[j])

实现非常直白：

lagrange_scale = np.empty_like(coords) for i in range(len(coords)): deltas = coords[i] - coords deltas[i] = 1.0 lagrange_scale[i] = 1.0 / np.prod(deltas)

为什么需要它？标准的 Lagrange 基函数

L_i(x) = ∏_{j≠i} (x − x_j) / (x_i − x_j)

其分母正是lagrange_scales返回值的倒数。把分母预计算成常量数组后，基函数在任意 x 处的求值就退化为一次多项式连乘再乘一个标量——这是 Warp 在 GPU 内核里高效求值形函数的惯用手段。

它的两个实际消费场景（源码证据）：

1. 四边形/六面体张量积形函数：SquareBipolynomialShapeFunctions 与 CubeTripolynomialShapeFunctions 构造时，先用quadrature_1d(point_count=degree+1, family=family)拿到节点（源码中变量名即lobatto_coords），再调lagrange_scales(lobatto_coords)得到LAGRANGE_SCALE，连同节点坐标、权重一起注册为wp.constant，供设备端基函数求值使用；
2. 节点坐标与节点权重：同一段代码把lobatto_coords/lobatto_weight存为LOBATTO_COORDS/LOBATTO_WEIGHT常量——注意这里节点与求积点是同一套点（求积点即插值节点），这正是 Gauss–Lobatto 节点家族的"双用途"设计：既精确积分，又天然适合做插值节点。

从一维到多维：张量积构造instantiate_quadrature

一维规则通过张量积扩展为二维/三维参考单元上的规则。PrototypeElement.instantiate_quadrature(order, family)（element.py）由Element各子类实现：

• LinearEdge：1D 规则原样返回（coords 补零成三元组）；
• Square：对 x、y 方向做张量积，weights = [wx * wy ...]，得到n×n个点；
• Cube：三方向张量积，weights = [wx * wy * wz ...]，得到n³个点。

# 在 2×2 网格上构造一个 2 阶 Gauss–Legendre 规则（每单元 2×2=4 点） geo = fem.Grid2D(res=wp.vec2i(2)) domain = fem.Cells(geo) quadrature = fem.RegularQuadrature(domain, order=2, family=fem.Polynomial.GAUSS_LEGENDRE)

注意当family=None时，instantiate_quadrature默认回退到GAUSS_LEGENDRE（element.py），这也是make_element_shape_function在未指定族时默认 LGL（shape/init.py）之外的另一个默认行为。

求积框架中的落地：RegularQuadrature与兄弟类

warp.fem的求积体系围绕Quadrature基类（quadrature.py）组织，它定义了point_count/point_coords/point_weight/point_index/point_evaluation_index等设备端接口，以及从"求积点求值索引 → 所属单元"的映射表。polynomial.quadrature_1d的直接消费者是RegularQuadrature：

• 构造时以(element, order, family, scalar_type)为键查CachedFormula缓存，命中则直接复用；
• 缓存内容由element.prototype.instantiate_quadrature(order, family)生成（即走本文所述张量积路径）；
• 点与权重被转成wp.array存入设备端Arg结构。源码注释特别说明：求积点/权重曾以 Warp 常量的形式传递，但点数多时容易引发寄存器溢出（register spilling），故改为数组参数（quadrature.py）——这是理解该设计演进的关键细节；
• 每单元点数N以wp.constant固化，point_index = N * domain_element_index + qp_index给出全局线性索引。

同框架下还有不依赖polynomial模块的兄弟类：NodalQuadrature（以空间节点为求积点，quadrature.py）、ExplicitQuadrature（用户逐单元给出点与权重，quadrature.py）以及测试中出现的PicQuadrature（粒子/质点求积）。它们在fem.integrate、fem.interpolate中作为quadrature=参数被消费。

测试与验证：单形积分如何证明规则正确

warp/tests/fem/test_fem_quadrature.py 中的test_regular_quadrature给出了规则的验证方法论：

1. 单形精度检验：对四种Polynomial族、阶数 0–7，用对应规则积分单项式x^degree，与解析值1/(degree+1)对比（places=4）。Gauss 族在高阶上依然精确，等距族在阶数逼近其精度上限时开始偏差；
2. 三角形上的变换检验：把三角形映射到方形后积分y^k1 · (1−x)^k2，与解析解1/((k1+k2+2)(k2+1))对照，验证张量积规则在非张量积单元上的行为；
3. 集成验证：test_nodal_quadrature（2 阶 LGL 单元每单元 9 点，积分x³y³得1/16）、test_particle_quadratures（显式规则积分）等，覆盖了从 1D 规则到完整装配管线的链路。

这些测试同时也是理解各家族适用精度的最佳教材：被积函数阶数越高、对精度越敏感，越应选用 Gauss–Legendre；需要端点节点做连续装配时，选 Lobatto–Gauss–Legendre。

实际项目中的用法参考

Warp 自带 FEM 示例中大量使用fem.Polynomial指定求积族：

• example_burgers.py：DG 方法中用fem.Polynomial.LOBATTO_GAUSS_LEGENDRE构造order=3的求积——DG 需要单元边界通量，LGL 的端点节点在此场景尤其合适；
• example_convection_diffusion_dg.py：同样选择 LGL 族；
• example_mixed_elasticity.py 与 example_streamlines.py：选择GAUSS_LEGENDRE族。

小结

warp.fem.polynomial虽小，却是 Warp 有限元管线的"数学基座"：

• quadrature_1d(point_count, family)以[0,1]归一化形式提供 GL（1–5 点）、LGL（2–5 点）、等距闭型/开型 Newton–Cotes（2–8 / 1–7 点）四族一维求积规则，权重和为 1；
• lagrange_scales(coords)预计算 Lagrange 基函数分母，支撑方形/立方体张量积形函数的高效设备端求值；
• 一维规则经PrototypeElement.instantiate_quadrature张量积化为 2D/3D 规则，由RegularQuadrature缓存并注入设备端；
• 族的选择（is_closed）同时决定形函数节点的拓扑分类（顶点/边/内部）与求积精度，是理解 Warp FEM 空间离散的关键概念。

深入阅读建议：polynomial 源码、求积框架、参考单元实现、形函数构造、求积测试 以及 warp.fem 模块总览。

【免费下载链接】warpA Python framework for GPU-accelerated simulation, robotics, and machine learning.项目地址: https://gitcode.com/GitHub_Trending/warp/warp