三角形单元插值函数很多人的第一反应是那组背了又忘、忘了又背的公式一次单元是三个面积坐标二次单元是 $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$ 次完全多项式的维数也就是 $(n1)(n2)/2$。把每个节点的坐标代进去要求取值为 $\delta_{qr}$解一个线性方程组系数就出来了。这个方法一次单元用着很舒服二次单元开始就难受。问题出在方程组的条件数上。三角形的节点坐标是仿射缩放的节点越密范德蒙德型矩阵的行列式越小条件数按指数级恶化。二次单元6×6 的方程组用双精度还能勉强解出来三次单元10×10就开始出现有效数字丢失四次单元基本没法直接用浮点算。更麻烦的是解出来的系数没有物理含义你没法从系数形式判断哪里算错了只能靠代入节点验证出错之后排查成本极高。更深一层的问题是不通用。每换一个阶次、每换一种节点布置就得重新解一遍方程组代码没法复用。而划线法的通用公式一次写出来任意阶次、任意维度都是用同一行代码这一点在后面第 2.4 节会展开。1.3 划线法的几何直觉用三族直线围猎一个节点划线法的思路完全反过来不去解方程而是用几何办法把非目标节点全部排除掉。具体操作是这样的。把三角形的三条边各做 $n$ 等分过每个等分点画平行于另外两条边的直线三族各 $n1$ 条直线在三角形内部交织成一张网网线的交点就是该阶次单元的插值节点总数恰好 $(n1)(n2)/2$。这张网有个关键性质每一个节点都唯一地落在三族直线中各一条之上。现在盯住其中一个目标节点 $\mathbf{p}_q$。要构造它的形函数就得让这个函数在其余所有节点上都为零。而其余每一个节点都躺在这张三族直线的某一条上。于是思路就出来了找出若干条恰好穿过所有非目标节点、且不穿过目标节点的直线把它们的方程乘起来乘积自然在所有非目标节点处为零。最后再乘一个常数把目标节点处的值拉到 1。这就是划线法的全部逻辑。它把解线性方程组换成了找直线把代数问题变成了几何问题而几何问题的解是可以直接用眼睛看出来的。一次单元只需要看三条边二次单元需要看中位线三次单元需要看三分之一等分线规律非常稳固。注意这里说的直线在面积坐标下就是 $L_a c$ 形式的等值线。后面所有推导都在面积坐标框架内完成所以先把这套坐标系讲透。2. 划线法的核心机制拆解2.1 面积坐标是自带刻度的三向标尺面积坐标 $L_1,L_2,L_3$ 的定义很直接$L_a$ 等于点 $\mathbf{p}$ 与第 $a$ 个顶点对边围成的小三角形面积除以单元总面积。三个分量满足$$L_1L_2L_31,\qquad L_a\ge 0\ (\mathbf{p}\text{ 在单元内})$$它和直角坐标的关系是线性的$$x\sum_{a1}^{3}L_a x_a,\qquad y\sum_{a1}^{3}L_a y_a$$关键在于它的等值线形状。$L_1c$ 是一条平行于边 2-3 的直线$c$ 从 0 到 1 变化时这条直线从边 2-3 平移到顶点 1。所以面积坐标天然就是一把三向标尺三个分量分别对应三个方向上的到对边的相对距离。节点处的取值也极其干净。第 $a$ 个顶点处$L_a1$其余为零边 2-3 上的任意点$L_10$。把三角形做 $n$ 等分之后落在距顶点 1 为 $i/n$那条平行线上的所有点都满足 $L_1i/n$。也就是说如果用多重指标 $(i,j,k)$、$ijkn$ 来标记节点那么该节点的面积坐标就是$$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_1i/n$、$L_2j/n$、$L_3k/n$。要把其余节点全部排除就沿着三个方向各自找若干条平行线。沿 $L_1$ 方向看所有 $L_1i/n$ 且 $i\neq i$ 的节点都在别的平行线上。第一族要取的直线是$$L_10,\ \frac{1}{n},\ \frac{2}{n},\ \dots,\ \frac{i-1}{n}$$一共 $i$ 条。为什么到 $i-1$ 就停因为 $L_1i/n$ 的那些平行线上节点的 $i$ 都大于 $i$而它们的三重指标之和必须等于 $n$如果 $ii$那么 $jkjk$能不能保证这些节点一定被另外两族直线排除能但需要更细致的论证。实际构造时我们用的是一种降阶思路形函数在每条 $L_1m/n$$mi$的平行线上都为零这 $i$ 条线的乘积就是$$\prod_{m0}^{i-1}\left(L_1-\frac{m}{n}\right)\frac{1}{n^i}\prod_{m0}^{i-1}(nL_1-m)$$同理沿 $L_2$ 方向取 $j$ 条沿 $L_3$ 方向取 $k$ 条。三族乘起来得到一个次数为 $ijkn$ 的 $n$ 次齐次多项式。这个多项式的零集包含了一张完整的网网上恰好覆盖了除目标节点之外的所有节点——这一点可以用多重指标的字典序或者单纯形的组合性质证明实操中直接代入验证就够。这里有个非常漂亮的等价写法。定义广义二项式系数$$\binom{x}{m}\frac{x(x-1)(x-2)\cdots(x-m1)}{m!}$$它本身就是一个 $m$ 次多项式根恰好是 $x0,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_11$取值 $2\times 12$所以得除以 2最终 $N_1L_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!$ 分母正好和在 $xm$ 处取值为 1这个要求匹配。换个角度看$\binom{nL_1}{i}$ 在 $L_1i/n$ 处的值始终为 1与 $n$ 和 $i$ 都无关所以三个因子相乘在目标节点处必然是 1。实操中这个差别很关键。我见过不少手写代码在归一化这一步漏除或者除错症状是刚度矩阵量级整体偏移但结构看起来是对的因为每个单元都错得一样装配之后位移场整体缩放了一个常数。这种错误靠肉眼查矩阵根本查不出来只能靠单位分解性检测。如果一开始就用二项式系数形式写这个坑直接绕过。2.4 通用公式一行代码覆盖所有阶次和维度把上面的论证收拢得到划线法的通用表达式。设单元阶次为 $n$节点用多重指标 $(i,j,k)$ 标记$ijkn$则$$N_{ijk}\binom{nL_1}{i}\binom{nL_2}{j}\binom{nL_3}{k}$$这个式子的适用性远超三角形。四面体上把面积坐标换成体积坐标 $L_1,\dots,L_4$多重指标换成四元组 $ijkln$公式完全一样只是因子多一个。二维三角形、三维四面体、甚至一维线段退化成两个因子就是普通的拉格朗日插值基函数全都套得进这一个模子。这也是划线法比待定系数法强的地方——待定系数法每换一个维度都得重解一遍方程组通用公式写出来之后连分支都不用加。单位分解性在这个形式下有一个极简的证明。对任意实数 $a,b,c$有范德蒙德卷积$$\sum_{ijkn}\binom{a}{i}\binom{b}{j}\binom{c}{k}\binom{abc}{n}$$把 $anL_1$、$bnL_2$、$cnL_3$ 代进去右边是 $\binom{n(L_1L_2L_3)}{n}\binom{n}{n}1$。所以 $\sum N_{ijk}\equiv 1$ 自动成立不需要任何额外条件。这条恒等式在组合数学里是很基础的东西搬到这里来单位分解性就从需要验证的性质变成了构造本身就带的属性。线性完备性也能顺手验证。令 $n1$$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$ 方向取 $i1$ 条线$L_10$也就是边 2-3$L_2$ 方向取 $j0$ 条空乘$L_3$ 方向取 $k0$ 条。所以形函数只由一条直线决定$$N_1\binom{L_1}{1}L_1$$验证一下顶点 1 处 $L_11$取值为 1顶点 2、3 处 $L_10$取值为 0边 2-3 上处处为零。三条要求一次满足。几何上的说法更直白顶点 1 的形函数零线就是它的对边。三角形三个顶点各自的零线就是各自的对边三条边把三个节点的势力范围划分得干干净净。这也解释了一次单元为什么在公共边上退化成一维线性插值——在边 1-2 上$L_30$$N_1L_1$、$N_2L_2$、$N_30$正好是这条边上两个端点的一维线性形函数。这一阶的实操没有任何坑但有一个基础设定值得强调节点编号必须逆时针。如果顶点顺序反了面积坐标的行列式会变号后面所有雅可比都是负的刚度矩阵整体带负号求解器要么报错要么给出荒谬结果。这个坑在导入外部网格文件时特别常见。3.2 二次单元顶点与边中点的分工二次单元是划线法真正开始显威力的地方。三边各二等分过等分点画平行线节点总数 6 个三个顶点、三个边中点。用多重指标标记$n2$多重指标 $(i,j,k)$面积坐标几何位置形函数(2,0,0)(1,0,0)顶点 1$L_1(2L_1-1)$(1,1,0)(1/2,1/2,0)边 1-2 中点$4L_1L_2$(1,0,1)(1/2,0,1/2)边 1-3 中点$4L_1L_3$(0,2,0)(0,1,0)顶点 2$L_2(2L_2-1)$(0,1,1)(0,1/2,1/2)边 2-3 中点$4L_2L_3$(0,0,2)(0,0,1)顶点 3$L_3(2L_3-1)$先看顶点 1$i2$。$L_1$ 方向要取两条线$L_10$ 和 $L_11/2$。前者是边 2-3后者是连接边 1-2 中点和边 1-3 中点的中位线。这两条平行线一夹除了顶点 1 自身其余 5 个节点全部落在这两条线上顶点 2、顶点 3、边 2-3 中点在 $L_10$ 上边 1-2 中点、边 1-3 中点在 $L_11/2$ 上。$jk0$另外两族不用取线。所以$$N_1\binom{2L_1}{2}\frac{2L_1(2L_1-1)}{2}L_1(2L_1-1)$$注意一个反直觉的性质这个形函数在 $L_11/2$中位线处取 0在 $L_11/4$ 处取 $-1/8$。也就是说它在单元内部是有负值的。这直接影响后面质量矩阵的性质第 5.2 节细说。再看边 1-2 中点$(i,j,k)(1,1,0)$。$L_1$ 方向取 $i1$ 条线$L_10$$L_2$ 方向取 $j1$ 条线$L_20$$L_3$ 方向不取。两条线分别是边 2-3 和边 1-3它们的交点是顶点 3。检查其余节点顶点 1 在 $L_20$ 上顶点 2 在 $L_10$ 上顶点 3 在两条线的交点上边 1-3 中点在 $L_20$ 上边 2-3 中点在 $L_10$ 上。除目标节点外全部覆盖。所以$$N\binom{2L_1}{1}\binom{2L_2}{1}4L_1L_2$$在目标节点处 $L_1L_21/2$$4\times 1/41$归一化自动正确。几何上这个选择很有讲究顶点节点只需要一族平行线就够了两条因为它离对边最远边中点节点需要两族各一条因为它