划线法构造三角形单元插值函数：从面积坐标到形函数验证

“划线法构造三角形单元插值函数”这个说法，第一次听容易让人以为是画网格的技巧。它确实是在画线，但画的不是网格线，而是面积坐标的等值线——把三角形按节点位置切成几块，插值函数的“零线”自然浮出来，形函数几乎不用解方程就能写出来。我在刚接手一个平面应力算例时，带我的工程师没让我背 N1 = L1、N4 = 4L1L2 这些公式，而是丢给我一张纸：“先把这个三角形切四块，切完你就知道为什么了。”那个下午我画了十几张三角形，后来这套笨办法在六节点、十节点单元上救了我很多次。这篇东西写给两类人：一类是刚学有限元、被形函数公式表劝退的同学，另一类是自己写程序、需要手推或验证三角形单元插值函数的开发者。读完之后你应该能做到：拿到任意三角形节点布置，先用划线法定性判断形函数长什么样，再用面积坐标把它精确写出来，最后用几行代码验证对不对。

1. 划线法到底在划什么：先把几何直觉立起来

1.1 从一维线段说起：二次插值的“划线”原型

划线法的核心思想在一维线段上就能看得很清楚。假设一条线段上布置三个节点：左端 ξ=0、中点 ξ=1/2、右端 ξ=1，要做二次拉格朗日插值。对应的三个基函数是 N0 = 2(ξ-1/2)(ξ-1)、N1 = -4ξ(ξ-1)、N2 = 2ξ(ξ-1/2)。注意看 N0 的零点：ξ=1/2 和 ξ=1，恰好是另外两个节点的位置。也就是说，每个基函数的零点集合，就是“除自己之外所有同类节点的等值线”。在一维里，这些“线”退化成点；到了三角形上，它们变成一族平行于单元边的直线，这就是“划线”两个字的来历。

我把这个过程理解成“用节点当钉子，把函数的零点位置钉死”。你不需要先假设一个多项式再解方程组，而是反过来：先确定这个形函数必须在哪些位置为零，再确定它必须在哪个位置为一，中间的形状由零点重数决定。一维二次基在 ξ=0 处取 1，在其他两个节点处取 0，所以它必然是 (ξ-1/2)(ξ-1) 乘一个常数，代入 ξ=0 定出常数 2。整个过程没有任何联立方程，全是几何定位加一个标定。三角形单元上的划线法，就是把这套动作从一条线搬到三个方向上。

1.2 三角形面积坐标：划线法天然的坐标系

要在三角形上划线，得先有一套描述“线”的语言，面积坐标就是这套语言。取三角形三个顶点为节点 1、2、3，内部任一点 P 与节点 2、3 围成一个小三角形，面积记作 A1；类似地定义 A2、A3。总面积 A = A1 + A2 + A3，于是定义 L1 = A1/A、L2 = A2/A、L3 = A3/A，三者之和恒等于 1。这套坐标有个极其好用的性质：L1 在节点 1 处等于 1，在节点 2、3 处等于 0；L2、L3 同理。换句话说，面积坐标本身就是三个“天然的线性形函数”。

更关键的是等值线的形状。L1 = 常数的点集，是一条平行于边 23 的直线。L1 = 0 就是边 23 本身，L1 = 1 退化成节点 1 一个点，L1 = 1/2 就是平行于边 23 的中位线。所以“在三角形里划线”，划的就是 L1 = 1/2、L1 = 1/3、L2 = 1/3 这一族平行线。这些线把大三角形切成若干小三角形和四边形，每个节点落在哪块区域里，形函数的非零支撑就在哪块区域里。这就是划线法能“看”出形函数的原因：零点线的位置决定了形函数的支撑范围，支撑范围又决定了形函数的乘积形式。

1.3 划线法与待定系数法的取舍：什么时候用哪个

构造三角形单元插值函数，教科书上通常给两条路。一条是待定系数法：设 N1 = a + bx + cy + dxy + …，把各节点的坐标代进去，令 N1 在节点 1 等于 1、在其他节点等于 0，解一个线性方程组定系数。另一条就是划线法：直接利用面积坐标的零点性质，把形函数写成乘积形式，再用节点值标定系数。

3 节点线性单元用待定系数法还算轻松，3×3 方程组手算几分钟。到了 6 节点二次单元，未知数变成 6 个，方程组里出现 xy、x²、y² 项，手算就开始出错；10 节点三次单元是 10×10 方程组，系数矩阵条件数还不小，纯手算基本不现实。划线法的优势在这里非常明显：6 节点单元的六个形函数，熟练之后十分钟能全部写出来，而且每一项都自带“在别的节点为零”的结构，不需要事后验证。

但划线法不是万能的。它要求节点布置规则——角点、边中点、边三等分点、重心，位置必须和面积坐标的等分线对得上。工程中遇到过渡单元、畸变单元、或者自己定义的特殊节点分布，划线法就不好使了，还是得回到待定系数法或者程序自动生成。我的习惯是：规则单元手推用划线法，程序里用统一的多项式拟合或拉格朗日插值公式兜底，两套结果互相印证，哪边出错一眼就能看出来。

2. 面积坐标的计算与三角形几何量落地

2.1 三种等价定义与行列式写法

面积坐标有三种常见定义，用哪种取决于你手头有什么。第一种是几何定义，就是上面说的面积比，直观但手算时要算三次三角形面积。第二种是线性插值定义：L1 是 x、y 的线性函数，满足 L1(x1,y1)=1、L1(x2,y2)=0、L1(x3,y3)=0，三个条件定一个线性函数，正好唯一。第三种是重心坐标定义，把点 P 写成三个顶点的凸组合 P = L1·P1 + L2·P2 + L3·P3，且 L1+L2+L3=1。三种定义在数学上完全等价，工程实现里最稳妥的是行列式写法。

设三角形顶点按逆时针排列，构造矩阵 M = [[x1,y1,1],[x2,y2,1],[x3,y3,1]]，则 2A = det(M)。把 P 点坐标替换第 i 行，得到的行列式除以 det(M)，就是 L_i。这种写法有个好处：符号自动正确。节点逆时针排列时 det(M) 为正，L_i 在三角形内部为正；节点顺时针排列时 det(M) 为负，L_i 全部变号，你能立刻从结果里发现编号顺序错了。很多手写代码用显式公式算面积坐标，公式里的正负号靠记忆，换个节点顺序就悄悄出错，行列式写法把这个问题堵死了。

2.2 一个具体三角形的完整手算

拿一个最简单的三角形练手：节点 1 在 (0,0)，节点 2 在 (2,0)，节点 3 在 (0,1)。逆时针排列，面积 A = 1，2A = 2。按行列式算：

L1 对应把 M 的第一行换成 (x,y,1)，得到 det[[x,y,1],[2,0,1],[0,1,1]] = x(0-1) - y(2-0) + 1(2-0) = -x - 2y + 2，除以 2 得 L1 = 1 - x/2 - y。L2 把第二行替换，det[[0,0,1],[x,y,1],[0,1,1]] = x，除以 2 得 L2 = x/2。L3 = y。验算：节点 1 (0,0) 处 L1=1、L2=L3=0；节点 2 (2,0) 处 L2=1；节点 3 (0,1) 处 L3=1，全部对得上。

接着往二次单元走。6 节点单元的边中点坐标是：节点 4 在边 12 中点 (1,0)，节点 5 在边 23 中点 (1,0.5)，节点 6 在边 31 中点 (0,0.5)。用形函数 N4 = 4L1L2 代入节点 4：L1 = 1 - 0.5 - 0 = 0.5，L2 = 0.5，N4 = 4×0.5×0.5 = 1。代入节点 1：L2=0，N4=0。代入节点 2：L1=0，N4=0。代入节点 5：L1 = 1-0.5-0.5 = 0，N4=0。六个节点全部验一遍，只有节点 4 处为 1，其余为 0。这就是划线法写出来的形函数该有的样子，不需要解任何方程。

2.3 雅可比矩阵与面积符号的坑

从面积坐标回到物理坐标，靠的是等参变换 x = Σ N_i x_i、y = Σ N_i y_i。求导链式法则要求雅可比矩阵 J = [[∂x/∂L1, ∂y/∂L1],[∂x/∂L2, ∂y/∂L2]]。对 3 节点线性单元，N_i 就是 L_i，J 是常数矩阵，det(J) 恒等于 2A，验算起来很轻松。对 6 节点、10 节点单元，N_i 含 L 的二次、三次项，J 随积分点变化，每个高斯点都要重新算，这是高阶单元计算量大的直接原因。

这里有个特别容易踩的坑：det(J) 的符号必须和单元面积符号一致。节点逆时针排列，det(J) = 2A > 0，一切正常；节点顺时针排列，det(J) < 0，数值积分出来的刚度矩阵整体负号，组装之后总刚出现负对角线，求解器要么报错要么给出物理上荒谬的位移场。更隐蔽的情况是部分单元顺时针、部分逆时针，总刚既不正定也不负定，残差曲线来回震荡，查起来非常痛苦。我的做法是在读入网格后加一道预处理：逐单元算 det(J)，为负就交换两个节点并记一条警告日志。这道检查花不了几毫秒，能省掉后面几小时的排查。

提示：面积坐标的等值线是平行于对边的直线，这个性质只在平面三角形上严格成立。到了曲边等参单元，物理空间里的“等 L1 线”会弯曲，但面积坐标作为参考坐标仍然成立，形函数形式不变，变的只是雅可比矩阵。

3. 用划线法构造3节点线性三角形单元

3.1 三条边就是三条零线

3 节点线性单元是划线法最朴素的展示。节点 1 的形函数必须在节点 2、3 处为零，而节点 2、3 连成的就是边 23，这条边恰好是 L1 = 0 的等值线。所以 N1 必然含因子 L1。节点 1 处 L1 = 1，要让 N1 = 1，系数只能是 1。于是 N1 = L1，同理 N2 = L2、N3 = L3。整个过程一句话：对边就是零线，形函数就是对应的面积坐标。

这个结论简单到容易被轻视，但它解释了一个重要现象：线性三角形单元能精确表示线性场。因为 L1、L2、L3 都是 x、y 的线性函数，任何线性场 u = a + bx + cy 都可以写成 u = u1·L1 + u2·L2 + u3·L3，节点值一插就准。这也是分片检验（patch test）能通过的根本原因。反过来，如果单元里出现了二次项，线性三角形就表示不了，必须上 6 节点单元。

3.2 三项性质验证：单位分解、节点性、完备性

形函数写完之后必须验三件事，这三件事在后续所有高阶单元上都要重复做。第一是单位分解：N1+N2+N3 = L1+L2+L3 = 1，恒成立，说明刚体平移能被精确表示。第二是节点性：N_i 在节点 j 处等于 δ_ij，这是插值函数的基本要求，也是组装时“节点值直接就是系数”的前提。第三是完备性：Σ N_i x_i = x、Σ N_i y_i = y，说明几何映射和场插值用的是同一套基函数，等参元的一致性靠它保证。

三项性质里，单位分解最容易被忽略，也最容易在高阶单元上出错。6 节点单元的形函数加起来必须恒等于 1，如果不等于 1，说明某个形函数写错了系数，或者漏了一项。我一般在代码里随机撒 100 个内部点，算 ΣN_i，最大值和最小值与 1 的偏差超过 1e-12 就报警。这个检查成本极低，但能抓住绝大多数手推公式的低级错误。

3.3 线性单元的刚度矩阵组装流程

线性单元的应变矩阵 B 是常数矩阵，因为 N_i 对 x、y 的偏导都是常数。平面应力问题里 B 是 3×6 矩阵，D 是 3×3 弹性矩阵，刚度矩阵 K = BᵀDB·t·A，t 是厚度，A 是单元面积，一次矩阵乘法就完事，不需要数值积分。这也是线性三角形单元在工程上仍然大量使用的原因：快，而且结果对网格畸变不敏感。

但快是有代价的。线性三角形单元的位移场在单元内是线性的，应变是常数，相邻单元之间应变不连续，应力云图看起来一块一块的。要缓解这个问题，要么加密网格，要么换 6 节点单元。我的经验是：做趋势分析、初步设计、参数扫描，线性三角形足够；要看局部应力集中、做疲劳评估，老老实实上二次单元。划线法在这两种单元上都能用，只是 6 节点单元要多切几刀。

4. 划线法构造6节点二次三角形单元

4.1 先画中位线：把三角形切成四块

6 节点单元的节点布置是：三个角点 1、2、3，三条边的中点 4（边 12）、5（边 23）、6（边 31）。划线法的第一步是画三条中位线：L1 = 1/2、L2 = 1/2、L3 = 1/2。三条线在重心附近交出一个中央小三角形，加上三个角上的小三角形，整个大三角形被切成四块。含顶点 1 的那块区域满足 L1 ≥ 1/2，含顶点 2 的区域满足 L2 ≥ 1/2，含顶点 3 的区域满足 L3 ≥ 1/2，中央那块三个坐标都小于 1/2。

这个切分直接把角点形函数的支撑范围限制住了。N1 应该在 L1 = 1/2 这条线上为零（因为这条线穿过节点 4、6，这两个节点处 N1 必须为 0），同时在对边 L1 = 0 上为零（节点 2、3、5 都在上面）。两个零点条件意味着 N1 含因子 L1 和 (L1 - 1/2)，写成 N1 = c·L1·(L1 - 1/2)。代入节点 1 处 L1 = 1，要求 N1 = 1，得 c = 2，所以 N1 = L1(2L1 - 1)。同理 N2 = L2(2L2 - 1)、N3 = L3(2L3 - 1)。

4.2 角点节点与边中点节点的划线差异

边中点节点的划线逻辑和角点不一样，这是很多人第一次学的时候卡住的地方。以节点 4（边 12 中点）为例，它必须在节点 1、2、3、5、6 处为零。节点 1、2 都在边 12 上，但 N4 在节点 4 处要等于 1，所以它不能在整条边 12 上为零。真正能用的零线是另外两条边所在的直线：L1 = 0（边 23）和 L2 = 0（边 31）。节点 3、5 在 L1 = 0 上，节点 6、3 在 L2 = 0 上，节点 1 处 L2 = 0，节点 2 处 L1 = 0，全部覆盖。所以 N4 含因子 L1·L2，写成 N4 = c·L1·L2，代入节点 4 处 L1 = L2 = 1/2，得 c = 4，于是 N4 = 4L1L2。

对比一下：角点形函数用了“对边 + 中位线”两条零线，支撑被压缩在角部小三角形里；边中点形函数只用了“两条邻边”两条零线，支撑铺满整个单元内部。这个差异在图上非常直观——角点形函数像一个小鼓包挤在角上，边中点形函数像一座横跨整个单元的山丘，峰值在边中点。这也解释了为什么 6 节点单元的边上位移是二次曲线：边 12 上 L3 = 0，N1 = L1(2L1-1)、N2 = L2(2L2-1)、N4 = 4L1L2，其余三个形函数为零，正好是一维二次拉格朗日基，边中点节点自然成为曲线的控制点。

注意：N4 = 4L1L2 里的系数 4 不是随便来的，它由“节点 4 处 L1 = L2 = 1/2、要求 N4 = 1”唯一确定。换成边三等分点，系数会变成 9/2，后面讲三次单元时会看到。

4.3 10节点三次单元的划线延展

三次单元（10 节点）的划线更密：每个方向画两条三等分线，L1 = 1/3、2/3，三个方向共六条线，把三角形切成九块小三角形，加上三个角点、六条边上的内点、一个重心点，正好 10 个节点。角点形函数要在 L1 = 0、1/3、2/3 三条线上为零，写成 N1 = c·L1·(L1 - 1/3)·(L1 - 2/3)，代入 L1 = 1 定出 c = 9/2，整理得 N1 = (1/2)L1(3L1 - 1)(3L1 - 2)。

边内节点的形函数稍微绕一点。以边 12 上靠节点 1 的那个三分点（记作节点 4，位于 L1 = 2/3、L2 = 1/3）为例，它需要在 L1 = 0、L2 = 0 两条边上为零，还需要在边 12 上的另外两个节点（节点 1 和节点 5）处为零。边 12 上 L3 = 0，一维坐标是 L1，节点位置是 L1 = 1、2/3、1/3。要让函数在 L1 = 1 和 L1 = 1/3 处为零，因子是 (L1 - 1)(L1 - 1/3)，再乘上 L1、L2 保证在两条边上为零，最后用节点 4 的值定系数。整理之后得到 N4 = (9/2)L1L2(3L1 - 1)。节点 5 对称，N5 = (9/2)L1L2(3L2 - 1)。重心节点 10 的形函数要在三条边上都为零，因子是 L1L2L3，代入重心 L1 = L2 = L3 = 1/3 定出系数 27，得 N10 = 27L1L2L3。

三次单元的手推工作量明显上来了，但逻辑和二次完全一致：先数零线，再写因子，最后用节点值标定系数。真正需要小心的是节点编号顺序，边内两个三分点谁靠谁，编号错了形函数就全乱，这一点在 6.2 节细说。

4.4 高阶单元的协调性与病态风险

协调性是划线法自动带来的福利。6 节点单元相邻两条边共享边中点，边上形函数退化成同一组一维二次基，位移在公共边上连续；10 节点单元共享两个三分点，边上退化成同一组一维三次基，同样连续。这个性质不需要额外证明，看边上的形函数表达式就能确认：边上的形函数只由该边的端点和边内节点决定，第三个面积坐标在边上恒为零，不参与。这也是为什么划线法特别适合教学——它把协调性变成了一个肉眼可见的事实。

病态风险则要单独提醒。高阶单元的形函数在节点附近变化剧烈，质量矩阵的条件数随阶数上升。10 节点单元如果网格畸变严重，或者单元尺寸悬殊，质量矩阵求逆可能出现数值问题。我在一个热传导算例里用过 10 节点单元，网格里有一个长宽比接近 20 的细长单元，结果温度场在角点附近出现非物理振荡，换成 6 节点单元加局部加密后消失。经验是：高阶单元配高质量网格，畸变大的区域退回低阶单元，或者用四边形单元过渡，不要在畸变单元上硬上三次插值。

5. 编程实现与数值验证

5.1 Python实现面积坐标与形函数

下面这段代码是我平时用来验证形函数的模板，面积坐标用行列式算，符号自动正确，6 节点形函数按 4.2 节的公式写。跑一遍就能把节点性、单位分解、线性完备性全验完。

import numpy as np def barycentric(x, y, xy): """xy: 3x2 数组，节点逆时针排列""" M = np.hstack([xy, np.ones((3, 1))]) detM = np.linalg.det(M) # = 2A，逆时针为正 L = np.zeros(3) for i in range(3): j, k = (i + 1) % 3, (i + 2) % 3 Mi = np.array([[x, y, 1.0], M[j], M[k]]) L[i] = np.linalg.det(Mi) / detM return L def shape6(L): L1, L2, L3 = L return np.array([ L1 * (2 * L1 - 1), L2 * (2 * L2 - 1), L3 * (2 * L3 - 1), 4 * L1 * L2, 4 * L2 * L3, 4 * L3 * L1, ]) xy = np.array([[0.0, 0.0], [2.0, 0.0], [0.0, 1.0]]) nodes6 = np.array([[0,0],[2,0],[0,1],[1,0],[1,0.5],[0,0.5]]) Nmat = np.zeros((6, 6)) for a, (x, y) in enumerate(nodes6): Nmat[a, :] = shape6(barycentric(x, y, xy)) print(np.round(Nmat, 12))

输出应该是一个 6×6 的单位矩阵，非对角元在 1e-16 量级。如果哪一行不是单位向量，先检查节点编号顺序，再检查形函数系数。这个测试我建议做成单元测试，每次改形函数都跑一遍，比肉眼检查可靠得多。

5.2 五项数值验证清单

第一项是节点性，就是上面那个单位矩阵测试。第二项是单位分解，在单元内随机取点，算 ΣN_i 与 1 的偏差，正常情况应该在 1e-14 以内。第三项是线性完备性，用形函数插值 x 坐标和 y 坐标，结果应该精确等于该点的物理坐标，偏差同样在 1e-14 量级。第四项是边协调性，取边 12 上的点（L3 强制为 0），确认 N3、N5、N6 都为零，且 N1+N2+N4 = 1。第五项是分片检验，构造一个常应变场，用几个单元拼一块方形区域，看节点位移是否精确复现线性场。

这五项里，第三项最容易暴露形函数系数错误，第四项最容易暴露编号错误。我遇到过一个问题：N4 和 N6 的系数都写对了，但节点 4 和节点 6 的坐标在数组里放反了，节点性测试居然还能过（因为单位矩阵对行置换不敏感），是边协调性测试报的警。所以我一般把五项测试打包成一个函数，一次性跑完。

5.3 数值积分与刚度矩阵的落地

6 节点单元的刚度矩阵需要数值积分，三角形上的高斯积分点用面积坐标表示。常用的 3 点规则精度到二次，7 点规则精度到五次，对二次单元来说 3 点就够，但为了保险我一般用 7 点。积分公式是 ∫ f dA ≈ A·Σ w_i·f(L_i)，权重和为 1。

# 三角形 7 点高斯积分（精度到 5 次） gauss7 = np.array([ [1/3, 1/3, 1/3, 0.225], [0.05971587, 0.47014206, 0.47014206, 0.13239415], [0.47014206, 0.05971587, 0.47014206, 0.13239415], [0.47014206, 0.47014206, 0.05971587, 0.13239415], [0.79742699, 0.10128651, 0.10128651, 0.12593918], [0.10128651, 0.79742699, 0.10128651, 0.12593918], [0.10128651, 0.10128651, 0.79742699, 0.12593918], ])

权重加起来是 1，积分时乘上单元面积 A。形函数对 x、y 的偏导要先对 L 求导，再用雅可比逆矩阵变换。6 节点单元的雅可比不是常数，每个积分点单独算，代码量比线性单元多，但结构一样。这里有个解析公式可以交叉验证：∫_T L1^a L2^b L3^c dA = (a!b!c!)/(a+b+c+2)! × 2A。比如 ∫ L1²L2² dA = (2!2!0!)/(6!) × 2A = A/180，用这个数去对数值积分结果，能快速判断积分点选得够不够。

6. 常见问题与排查技巧实录

6.1 面积算成负数、形函数飞出天际

面积坐标算出来是负数，最常见的原因是节点顺序顺时针。行列式写法下 det(M) 为负，三个 L 全部变号，形函数在单元内不再单调，刚度矩阵出现负特征值。解决办法很简单：统一所有单元为逆时针编号，或者在读网格时检测 det(M) 的符号，为负就交换节点 2 和节点 3，同时记一条日志。不要图省事用绝对值，绝对值会让左右手系混在一起，应力符号会出错。

形函数数值“飞出天际”，比如在单元内部算出 N 的绝对值远大于 1，通常是三角形退化了——三个节点接近共线，面积接近零，面积坐标的分母极小，数值爆炸。这种情况要先查网格质量，算一下每个单元的面积和最小高，面积小于平均面积千分之一的单元单独标出来。我一般要求最小高不小于最长边的 1/20，低于这个值就局部加密。退化的三角形不仅形函数难看，雅可比也接近奇异，刚度矩阵条件数会失控。

还有一种情况是面积坐标之和明显偏离 1，比如 1.0001。这多半是浮点累积误差，量级在 1e-12 以内可以忽略；如果到了 1e-6，说明你的坐标数据精度不够，或者混用了单精度和双精度。有限元程序里所有几何量一律用双精度，单精度省不了多少内存，但能在这种地方坑你半天。

6.2 节点编号逆序导致的雅可比翻负

6 节点单元的节点编号有个固定约定：1、2、3 是角点，4 是边 12 中点，5 是边 23 中点，6 是边 31 中点。这个顺序一旦打乱，形函数和节点坐标的对应关系就错位，程序不会报错，但结果会悄悄偏掉。我见过一个算例，节点 4 和节点 6 在连接表里写反了，线性场测试还能勉强通过（误差 1e-3 量级），二次场测试直接崩掉，位移场在边上出现折点。排查方法就是 5.2 节的边协调性测试：取边 12 中点，确认 N4 = 1、其余为 0；如果 N6 = 1，说明编号反了。

10 节点单元更麻烦，每条边上有两个三分点，谁靠谁必须和形函数公式里的 L1、L2 对应。我的做法是在单元定义里显式存一个节点局部坐标表，组装前用形函数在局部坐标处求值，确认得到单位矩阵，再进入正式计算。这张表看起来冗余，但它是唯一能防住编号错位的硬约束。

提示：不同教材对 10 节点单元的编号约定不完全一样，有的是边 12 上先 4 后 5、边 23 上先 6 后 7，有的是按逆时针连续编号。抄公式之前先确认编号约定，不要公式抄一半、编号用另一半。

6.3 常见问题速查表

这张表是我这些年踩坑踩出来的，基本上形函数相关的问题都能在里面找到影子。真正需要记住的是两条：符号问题查编号方向，数值问题查网格质量。大部分“形函数不对”的报错，最后追下去都不是形函数公式写错，而是节点顺序或者坐标数据的问题。

最后分享一个我自己的习惯。每次手推完一组新的形函数，我不急着往程序里塞，而是先在纸上画三条线，标出每个节点的零线位置，然后用 5.1 节那段代码跑一遍单位矩阵测试。纸上的图和代码的结果对得上，才继续往下做。这个流程看起来慢，但比在几百行组装代码里找一个小数点要快得多。划线法的价值也在这里：它让你在动手算之前，先在脑子里把形函数的形状看见。