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_scalesLagrange 多项式的缩放因子quadrature_1d一维求积点与权重。模块本身是薄封装层warp/fem/polynomial.py 仅做重新导出真正的实现位于 warp/_src/fem/polynomial.pyfrom 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它决定了一维区间上的插值节点/求积节点布局。源码 定义了四种枚举成员字符串值是否含端点典型用途Polynomial.GAUSS_LEGENDREGL否精确度高、不含端点的经典高斯求积Polynomial.LOBATTO_GAUSS_LEGENDRELGL是含端点利于节点与单元边界对齐C⁰ 连续装配Polynomial.EQUISPACED_CLOSEDclosed是等距闭型 Newton–Cotes梯形、Simpson 等Polynomial.EQUISPACED_OPENopen否等距开型 Newton–Cotes辅助函数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一维求积规则 APIquadrature_1d(point_count, family)返回(coords, weights)二元组coordsNumPy 数组求积点坐标统一归一化到[0, 1]区间weightsNumPy 数组对应权重按区间长度归一化∫₀¹ 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硬编码了 n1 到 n5 的经典高斯节点内部先在[-1, 1]上构造再整体移位缩放至[0, 1]n节点[0,1] 区间权重10.51.020.21132487, 0.788675130.5, 0.530.5, 0.11270167, 0.887298334/9, 5/18, 5/1840.06943184, 0.33000948, 0.66999052, 0.93056816(18√30)/72, (18−√30)/72对称两对50.5 及两对对称节点128/450 及 (322±13√70)/1800对称两对n 个 Gauss–Legendre 点对不超过 2n−1 阶多项式精确成立n3 时权重和为 1可精确积分到 5 阶是精度/点数比最高的规则。它的代价是不含端点——用于装配时需要在单元边界做额外的通量处理例如 DG 方法。Lobatto–Gauss–Legendre 族_lobatto_gauss_legendre_quadrature_1d支持 n2 到 n5强制包含区间两端点n节点权重20, 10.5, 0.5即梯形法则30, 0.5, 11/6, 2/3, 1/6即 Simpson 法则40, 0.27639320, 0.72360680, 11/12, 5/12, 5/12, 1/1250, 0.17267316, 0.5, 0.82732684, 11/20, 49/180, 16/45, 49/180, 1/20LGL 的端点节点使其成为连续有限元如 C⁰ Lagrange的天然选择节点与单元顶点重合方便全局自由度编号与装配。其精度为 2n−3 阶n3 时对应 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。闭型代表规则归一化后n2[0.5, 0.5]梯形法则1 阶精度n3[1/6, 2/3, 1/6]Simpson 法则3 阶精度n4[1/8, 3/8, 3/8, 1/8]Simpson 3/8 法则n5[14/180, 64/180, 24/180, 64/180, 14/180]Boole 法则更高阶 n 直至 8 均有硬编码权重如 n8 时中心权重含343/640。开型族节点避开端点则出现负权重n3 时[2, −1, 2]/3、n5 时[11, −14, 26, −14, 11]/20。负权重意味着积分结果对舍入误差更敏感高次等距 Newton–Cotes 也因 Runge 现象而不适合高阶使用——这正是源码将高阶精度任务交给 Gauss 族的原因。需要指出的是等距族在数值上不如高斯族稳定实践中主要服务于等距网格与简单测试场景。支持的点数范围四种族的实现覆盖点数有限超出即抛NotImplementedErrorGL1–5 点LGL2–5 点闭型 Newton–Cotes2–8 点开型 Newton–Cotes1–7 点。从阶数到点数_point_count_from_order的换算规则quadrature_1d直接吃点数而有限元语境通常从多项式阶数出发。参考单元层在 element.py 中提供了_point_count_from_order(order, family)完成换算family点数公式说明GAUSS_LEGENDREmax(1, order // 2 1)每 2 阶加 1 个点LOBATTO_GAUSS_LEGENDREmax(2, order // 2 2)端点占 2 个点其余每 2 阶加 1EQUISPACED_CLOSEDmax(2, 2 * (order // 2) 1)闭型奇数个等距点EQUISPACED_OPENmax(1, 2 * (order // 2) 1)开型奇数个等距点例如order2、familyGAUSS_LEGENDRE得到 2 个点两点 Gauss 规则可精确积分 3 阶多项式满足二次被积函数需要而order2的 LGL 需要2//22 3个点Simpson 法则。这一映射保证了积分精度不低于被积多项式的阶数。lagrange_scalesLagrange 形函数的缩放因子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 内核里高效求值形函数的惯用手段。它的两个实际消费场景源码证据四边形/六面体张量积形函数SquareBipolynomialShapeFunctions 与 CubeTripolynomialShapeFunctions 构造时先用quadrature_1d(point_countdegree1, familyfamily)拿到节点源码中变量名即lobatto_coords再调lagrange_scales(lobatto_coords)得到LAGRANGE_SCALE连同节点坐标、权重一起注册为wp.constant供设备端基函数求值使用节点坐标与节点权重同一段代码把lobatto_coords/lobatto_weight存为LOBATTO_COORDS/LOBATTO_WEIGHT常量——注意这里节点与求积点是同一套点求积点即插值节点这正是 Gauss–Lobatto 节点家族的双用途设计既精确积分又天然适合做插值节点。从一维到多维张量积构造instantiate_quadrature一维规则通过张量积扩展为二维/三维参考单元上的规则。PrototypeElement.instantiate_quadrature(order, family)element.py由Element各子类实现LinearEdge1D 规则原样返回coords 补零成三元组Square对 x、y 方向做张量积weights [wx * wy ...]得到n×n个点Cube三方向张量积weights [wx * wy * wz ...]得到n³个点。# 在 2×2 网格上构造一个 2 阶 Gauss–Legendre 规则每单元 2×24 点 geo fem.Grid2D(reswp.vec2i(2)) domain fem.Cells(geo) quadrature fem.RegularQuadrature(domain, order2, familyfem.Polynomial.GAUSS_LEGENDRE)注意当familyNone时instantiate_quadrature默认回退到GAUSS_LEGENDREelement.py这也是make_element_shape_function在未指定族时默认 LGLshape/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给出了规则的验证方法论单形精度检验对四种Polynomial族、阶数 0–7用对应规则积分单项式x^degree与解析值1/(degree1)对比places4。Gauss 族在高阶上依然精确等距族在阶数逼近其精度上限时开始偏差三角形上的变换检验把三角形映射到方形后积分y^k1 · (1−x)^k2与解析解1/((k1k22)(k21))对照验证张量积规则在非张量积单元上的行为集成验证test_nodal_quadrature2 阶 LGL 单元每单元 9 点积分x³y³得1/16、test_particle_quadratures显式规则积分等覆盖了从 1D 规则到完整装配管线的链路。这些测试同时也是理解各家族适用精度的最佳教材被积函数阶数越高、对精度越敏感越应选用 Gauss–Legendre需要端点节点做连续装配时选 Lobatto–Gauss–Legendre。实际项目中的用法参考Warp 自带 FEM 示例中大量使用fem.Polynomial指定求积族example_burgers.pyDG 方法中用fem.Polynomial.LOBATTO_GAUSS_LEGENDRE构造order3的求积——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]归一化形式提供 GL1–5 点、LGL2–5 点、等距闭型/开型 Newton–Cotes2–8 / 1–7 点四族一维求积规则权重和为 1lagrange_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创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考