有限元三角形单元形函数划线法：通用公式与高阶单元实现

三角形单元插值函数，很多人的第一反应是那组背了又忘、忘了又背的公式：一次单元是三个面积坐标，二次单元是 $L_1(2L_1-1)$ 和 $4L_1L_2$，三次单元往后就是一串带 $9/2$、$27/2$ 系数的多项式，看着像天书。划线法解决的恰恰是这个痛点——它不要求你记住任何一条形函数，只要在三角形里把平行于三条边的等分线画出来，形函数自己就"浮"出来了。这套方法在有限元教材里常常只用半页纸带过，但它背后牵扯的东西不少：零集分离的构造逻辑、二项式系数形式的通用公式、范德蒙德卷积给出的单位分解证明，以及一整套高阶单元落地时踩过的坑。下面我按自己推公式、写代码、调程序的顺序，把这条线从头到尾捋一遍，从一次单元一路推到三次单元和四面体单元，每个关键参数都给出计算过程，代码可以直接跑。

1. 三角形单元插值函数到底要满足什么

1.1 三个硬条件：插值性、单位分解、分片连续

先说清楚目标函数长什么样。给定一个三角形单元，顶点记为 1、2、3，内部布了若干个插值节点 $\mathbf{p}_q$，我们要找一组多项式 $N_q(\mathbf{x})$，让它们替我们"搬运"节点上的物理量。这组多项式必须同时满足三条硬性要求，缺一条整个有限元框架就散架。

第一条是插值性，业内也叫 Kronecker delta 性质：$N_q(\mathbf{p}r)=\delta{qr}$。也就是第 $q$ 个形函数在自己节点上取 1，在别人节点上取 0。这条要求看起来朴素，实际是整个有限元后处理的根基——只有它成立，单元内的位移场才能写成 $\mathbf{u}=\sum_q N_q \mathbf{u}_q$，节点自由度才有明确的物理含义，边界条件的施加才能一个节点一个节点地"钉"。

第二条是单位分解性，$\sum_q N_q \equiv 1$。这条保证常应变状态能被精确表示：所有节点给同样的位移，单元内每一点的位移就该等于这个常数。如果形函数之和不是 1，单元在刚体平移下就会产生虚假应变，网格细化之后位移不收敛。

第三条是分片连续性，相邻单元在公共边上的插值结果必须一致。三角形单元的边是直线，公共边上只有该边上的节点共享，所以只要形函数在边上退化成一维插值函数，而且两侧用的是同一组边节点，连续性自动满足。

这三条里，第一条最容易满足（乘个系数就行），第二条最容易被忽略（高阶单元手推公式时经常漏掉归一化），第三条最能决定构造方法的选择——它逼着你必须用"边界退化为低维"的多项式，而三族平行直线恰好提供了这个性质。

1.2 待定系数法为什么在高阶单元上不好用

最朴素的构造思路是待定系数法：设形函数为完全 $n$ 次多项式，未知系数个数等于三角形上 $n$ 次完全多项式的维数，也就是 $(n+1)(n+2)/2$。把每个节点的坐标代进去，要求取值为 $\delta_{qr}$，解一个线性方程组，系数就出来了。

这个方法一次单元用着很舒服，二次单元开始就难受。问题出在方程组的条件数上。三角形的节点坐标是仿射缩放的，节点越密，范德蒙德型矩阵的行列式越小，条件数按指数级恶化。二次单元（6×6 的方程组）用双精度还能勉强解出来，三次单元（10×10）就开始出现有效数字丢失，四次单元基本没法直接用浮点算。更麻烦的是，解出来的系数没有物理含义，你没法从系数形式判断哪里算错了，只能靠代入节点验证，出错之后排查成本极高。

更深一层的问题是不通用。每换一个阶次、每换一种节点布置，就得重新解一遍方程组，代码没法复用。而划线法的通用公式一次写出来，任意阶次、任意维度都是用同一行代码，这一点在后面第 2.4 节会展开。

1.3 划线法的几何直觉：用三族直线"围猎"一个节点

划线法的思路完全反过来：不去解方程，而是用几何办法把"非目标节点"全部排除掉。

具体操作是这样的。把三角形的三条边各做 $n$ 等分，过每个等分点画平行于另外两条边的直线，三族各 $n+1$ 条直线在三角形内部交织成一张网，网线的交点就是该阶次单元的插值节点，总数恰好 $(n+1)(n+2)/2$。这张网有个关键性质：每一个节点，都唯一地落在三族直线中各一条之上。

现在盯住其中一个目标节点 $\mathbf{p}_q$。要构造它的形函数，就得让这个函数在其余所有节点上都为零。而其余每一个节点，都躺在这张三族直线的某一条上。于是思路就出来了：找出若干条"恰好穿过所有非目标节点、且不穿过目标节点"的直线，把它们的方程乘起来，乘积自然在所有非目标节点处为零。最后再乘一个常数，把目标节点处的值拉到 1。

这就是划线法的全部逻辑。它把"解线性方程组"换成了"找直线"，把代数问题变成了几何问题，而几何问题的解是可以直接用眼睛看出来的。一次单元只需要看三条边，二次单元需要看中位线，三次单元需要看三分之一等分线，规律非常稳固。

注意：这里说的"直线"，在面积坐标下就是 $L_a = c$ 形式的等值线。后面所有推导都在面积坐标框架内完成，所以先把这套坐标系讲透。

2. 划线法的核心机制拆解

2.1 面积坐标是自带刻度的三向标尺

面积坐标 $L_1,L_2,L_3$ 的定义很直接：$L_a$ 等于点 $\mathbf{p}$ 与第 $a$ 个顶点对边围成的小三角形面积，除以单元总面积。三个分量满足

$$L_1+L_2+L_3=1,\qquad L_a\ge 0\ (\mathbf{p}\text{ 在单元内})$$

它和直角坐标的关系是线性的：

$$x=\sum_{a=1}^{3}L_a x_a,\qquad y=\sum_{a=1}^{3}L_a y_a$$

关键在于它的等值线形状。$L_1=c$ 是一条平行于边 2-3 的直线，$c$ 从 0 到 1 变化时，这条直线从边 2-3 平移到顶点 1。所以面积坐标天然就是一把三向标尺：三个分量分别对应三个方向上的"到对边的相对距离"。

节点处的取值也极其干净。第 $a$ 个顶点处，$L_a=1$，其余为零；边 2-3 上的任意点，$L_1=0$。把三角形做 $n$ 等分之后，落在"距顶点 1 为 $i/n$"那条平行线上的所有点，都满足 $L_1=i/n$。也就是说，如果用多重指标 $(i,j,k)$、$i+j+k=n$ 来标记节点，那么该节点的面积坐标就是

$$L_1=\frac{i}{n},\qquad L_2=\frac{j}{n},\qquad L_3=\frac{k}{n}$$

这套记号是整个划线法能写成通用公式的前提。看到 $(i,j,k)$，立刻知道它在三族直线网里的位置，不用画图也能推。

2.2 零集分离：形函数的骨架由直线族决定

现在把"排除非目标节点"这件事写严谨。目标节点是 $(i,j,k)$，也就是 $L_1=i/n$、$L_2=j/n$、$L_3=k/n$。要把其余节点全部排除，就沿着三个方向各自找若干条平行线。

沿 $L_1$ 方向看：所有 $L_1=i'/n$ 且 $i'\neq i$ 的节点都在别的平行线上。第一族要取的直线是

$$L_1=0,\ \frac{1}{n},\ \frac{2}{n},\ \dots,\ \frac{i-1}{n}$$

一共 $i$ 条。为什么到 $i-1$ 就停？因为 $L_1>i/n$ 的那些平行线上，节点的 $i'$ 都大于 $i$，而它们的三重指标之和必须等于 $n$，如果 $i'>i$，那么 $j'+k'<j+k$，能不能保证这些节点一定被另外两族直线排除？能，但需要更细致的论证。实际构造时我们用的是一种"降阶"思路：形函数在每条 $L_1=m/n$（$m<i$）的平行线上都为零，这 $i$ 条线的乘积就是

$$\prod_{m=0}^{i-1}\left(L_1-\frac{m}{n}\right)=\frac{1}{n^i}\prod_{m=0}^{i-1}(nL_1-m)$$

同理，沿 $L_2$ 方向取 $j$ 条，沿 $L_3$ 方向取 $k$ 条。三族乘起来，得到一个次数为 $i+j+k=n$ 的 $n$ 次齐次多项式。这个多项式的零集包含了一张完整的"网"，网上恰好覆盖了除目标节点之外的所有节点——这一点可以用多重指标的字典序或者单纯形的组合性质证明，实操中直接代入验证就够。

这里有个非常漂亮的等价写法。定义广义二项式系数

$$\binom{x}{m}=\frac{x(x-1)(x-2)\cdots(x-m+1)}{m!}$$

它本身就是一个 $m$ 次多项式，根恰好是 $x=0,1,\dots,m-1$。把上面那串乘积除以 $i!$，正好凑成 $\binom{nL_1}{i}$。于是零集多项式可以直接写成

$$\Phi_{ijk}=\binom{nL_1}{i}\binom{nL_2}{j}\binom{nL_3}{k}$$

这个形式的妙处马上就能看到：在目标节点 $(i,j,k)$ 处，$\binom{nL_1}{i}=\binom{i}{i}=1$，同理另外两个也等于 1，所以 $\Phi_{ijk}$ 在目标节点处的取值恰好就是 1。也就是说，归一化常数自动就是 1，一个多余的除法都不用做。这一点如果用原始的乘积形式写，就需要显式地算一遍 $\prod$ 在节点处的值再相除，代码里还容易因为浮点误差把 1 除成 1.0000000002。

2.3 归一化系数为什么能"天然等于1"

上一节结论值得单独展开，因为它是划线法区别于"凑公式"的关键。

用原始乘积形式，构造出来的是 $\tilde N_{ijk}=\prod(nL_1-m)\prod(nL_2-m)\prod(nL_3-m)$，这个式子只保证在非目标节点为零，在目标节点处取什么值完全看运气。以二次单元的顶点 1 为例，$(i,j,k)=(2,0,0)$，得到

$$\tilde N= (2L_1)(2L_1-1)$$

在顶点 1 处 $L_1=1$，取值 $2\times 1=2$，所以得除以 2，最终 $N_1=L_1(2L_1-1)$。这就是那个"莫名其妙多出来的系数"的来源——它其实是归一化因子。

换成二项式系数形式，$\binom{2L_1}{2}=\frac{2L_1(2L_1-1)}{2}=L_1(2L_1-1)$，一步到位，系数 1。这不是巧合：$\binom{x}{m}$ 的定义里那个 $m!$ 分母，正好和"在 $x=m$ 处取值为 1"这个要求匹配。换个角度看，$\binom{nL_1}{i}$ 在 $L_1=i/n$ 处的值始终为 1，与 $n$ 和 $i$ 都无关，所以三个因子相乘在目标节点处必然是 1。

实操中这个差别很关键。我见过不少手写代码在归一化这一步漏除或者除错，症状是刚度矩阵量级整体偏移，但结构看起来"是对的"，因为每个单元都错得一样，装配之后位移场整体缩放了一个常数。这种错误靠肉眼查矩阵根本查不出来，只能靠单位分解性检测。如果一开始就用二项式系数形式写，这个坑直接绕过。

2.4 通用公式：一行代码覆盖所有阶次和维度

把上面的论证收拢，得到划线法的通用表达式。设单元阶次为 $n$，节点用多重指标 $(i,j,k)$ 标记，$i+j+k=n$，则

$$N_{ijk}=\binom{nL_1}{i}\binom{nL_2}{j}\binom{nL_3}{k}$$

这个式子的适用性远超三角形。四面体上把面积坐标换成体积坐标 $L_1,\dots,L_4$，多重指标换成四元组 $i+j+k+l=n$，公式完全一样，只是因子多一个。二维三角形、三维四面体、甚至一维线段（退化成两个因子，就是普通的拉格朗日插值基函数），全都套得进这一个模子。这也是划线法比待定系数法强的地方——待定系数法每换一个维度都得重解一遍方程组，通用公式写出来之后连分支都不用加。

单位分解性在这个形式下有一个极简的证明。对任意实数 $a,b,c$，有范德蒙德卷积

$$\sum_{i+j+k=n}\binom{a}{i}\binom{b}{j}\binom{c}{k}=\binom{a+b+c}{n}$$

把 $a=nL_1$、$b=nL_2$、$c=nL_3$ 代进去，右边是 $\binom{n(L_1+L_2+L_3)}{n}=\binom{n}{n}=1$。所以 $\sum N_{ijk}\equiv 1$ 自动成立，不需要任何额外条件。这条恒等式在组合数学里是很基础的东西，搬到这里来，单位分解性就从"需要验证的性质"变成了"构造本身就带的属性"。

线性完备性也能顺手验证。令 $n=1$，$N_{100}=L_1$、$N_{010}=L_2$、$N_{001}=L_3$，而 $x=\sum L_a x_a$ 是三个形函数的线性组合，说明一次单元能精确表示任意线性场。高阶单元自然包含所有低阶多项式，完备性逐阶继承。

3. 从一次到三次的逐阶实操

3.1 一次单元：三条边线一画就完事

一阶是最简单的情形，但把流程完整走一遍有助于理解后面的推广。

三边不做等分，节点就是三个顶点。针对顶点 1，$(i,j,k)=(1,0,0)$。按 2.2 节的规则，$L_1$ 方向取 $i=1$ 条线：$L_1=0$，也就是边 2-3；$L_2$ 方向取 $j=0$ 条（空乘），$L_3$ 方向取 $k=0$ 条。所以形函数只由一条直线决定：

$$N_1=\binom{L_1}{1}=L_1$$

验证一下：顶点 1 处 $L_1=1$，取值为 1；顶点 2、3 处 $L_1=0$，取值为 0；边 2-3 上处处为零。三条要求一次满足。

几何上的说法更直白：顶点 1 的形函数，"零线"就是它的对边。三角形三个顶点各自的零线就是各自的对边，三条边把三个节点的"势力范围"划分得干干净净。这也解释了一次单元为什么在公共边上退化成一维线性插值——在边 1-2 上，$L_3=0$，$N_1=L_1$、$N_2=L_2$、$N_3=0$，正好是这条边上两个端点的一维线性形函数。

这一阶的实操没有任何坑，但有一个基础设定值得强调：节点编号必须逆时针。如果顶点顺序反了，面积坐标的行列式会变号，后面所有雅可比都是负的，刚度矩阵整体带负号，求解器要么报错要么给出荒谬结果。这个坑在导入外部网格文件时特别常见。

3.2 二次单元：顶点与边中点的分工

二次单元是划线法真正开始显威力的地方。三边各二等分，过等分点画平行线，节点总数 6 个：三个顶点、三个边中点。用多重指标标记，$n=2$：

先看顶点 1，$i=2$。$L_1$ 方向要取两条线：$L_1=0$ 和 $L_1=1/2$。前者是边 2-3，后者是连接边 1-2 中点和边 1-3 中点的中位线。这两条平行线一夹，除了顶点 1 自身，其余 5 个节点全部落在这两条线上：顶点 2、顶点 3、边 2-3 中点在 $L_1=0$ 上；边 1-2 中点、边 1-3 中点在 $L_1=1/2$ 上。$j=k=0$，另外两族不用取线。所以

$$N_1=\binom{2L_1}{2}=\frac{2L_1(2L_1-1)}{2}=L_1(2L_1-1)$$

注意一个反直觉的性质：这个形函数在 $L_1=1/2$（中位线）处取 0，在 $L_1=1/4$ 处取 $-1/8$。也就是说它在单元内部是有负值的。这直接影响后面质量矩阵的性质，第 5.2 节细说。

再看边 1-2 中点，$(i,j,k)=(1,1,0)$。$L_1$ 方向取 $i=1$ 条线：$L_1=0$；$L_2$ 方向取 $j=1$ 条线：$L_2=0$；$L_3$ 方向不取。两条线分别是边 2-3 和边 1-3，它们的交点是顶点 3。检查其余节点：顶点 1 在 $L_2=0$ 上，顶点 2 在 $L_1=0$ 上，顶点 3 在两条线的交点上，边 1-3 中点在 $L_2=0$ 上，边 2-3 中点在 $L_1=0$ 上。除目标节点外全部覆盖。所以

$$N=\binom{2L_1}{1}\binom{2L_2}{1}=4L_1L_2$$

在目标节点处 $L_1=L_2=1/2$，$4\times 1/4=1$，归一化自动正确。

几何上这个选择很有讲究：顶点节点只需要一族平行线就够了（两条），因为它离对边最远；边中点节点需要两族（各一条），因为它