这篇笔记在做什么

这一篇把原始对话里关于二维标量有限元的内容整理成一条完整链：控制方程、弱形式、全局编号、局部基函数、三角形单元矩阵、边界线段矩阵、正方形网格算例，以及和三维电磁棱边元的对应关系。

如果你想先在二维问题上把有限元的骨架走通，再回到三维电磁问题，这篇最适合当中间桥梁。

二维 Helmholtz 方程有限元方法

三角形单元 · 线性基函数 · Robin 边界条件 · 完整推导

一、问题陈述

1.1 控制方程

求 $\phi(x,y)$，满足：

$$ -\frac{\partial}{\partial x}\left(\alpha_x \frac{\partial \phi}{\partial x}\right)-\frac{\partial}{\partial y}\left(\alpha_y \frac{\partial \phi}{\partial y}\right)+\beta\phi=f, \quad (x,y)\in\Omega \tag{1} $$

1.2 Robin 边界条件

$$ \left(\alpha_x \frac{\partial \phi}{\partial x}\hat{x}+\alpha_y \frac{\partial \phi}{\partial y}\hat{y}\right)\cdot\hat{n}+\gamma \phi=q, \quad (x,y)\in\Gamma_2=\partial\Omega \tag{2} $$

1.3 记号简化

定义各向异性梯度算子：

$$ \bar{\bar{\alpha}}\nabla\phi \equiv \alpha_x \frac{\partial \phi}{\partial x}\hat{x}+\alpha_y \frac{\partial \phi}{\partial y}\hat{y} $$

其中 $\bar{\bar{\alpha}} = \text{diag}(\alpha_x, \alpha_y)$。

则方程和边界条件简写为：

$$ -\nabla\cdot(\bar{\bar{\alpha}}\nabla\phi)+\beta\phi = f \quad \text{in } \Omega \tag{1'} $$$$ (\bar{\bar{\alpha}}\nabla\phi)\cdot\hat{n}+\gamma\phi = q \quad \text{on } \Gamma_2 \tag{2'} $$

二、弱形式（变分形式）的推导

2.1 加权残量

用测试函数 $w(x,y)$ 乘以方程 (1)，在 $\Omega$ 上积分：

$$ \int_\Omega w\left[-\frac{\partial}{\partial x}\left(\alpha_x \frac{\partial \phi}{\partial x}\right)-\frac{\partial}{\partial y}\left(\alpha_y \frac{\partial \phi}{\partial y}\right)+\beta\phi\right]d\Omega = \int_\Omega w\,f\,d\Omega $$

简写为：

$$ -\int_\Omega w\,\nabla\cdot(\bar{\bar{\alpha}}\nabla\phi)\,d\Omega + \int_\Omega w\,\beta\,\phi\,d\Omega = \int_\Omega w\,f\,d\Omega \tag{3} $$

2.2 对第一项进行分部积分

利用矢量恒等式：

$$ w\,\nabla\cdot\mathbf{F} = \nabla\cdot(w\mathbf{F}) - \nabla w\cdot\mathbf{F} $$

取 $\mathbf{F} = \bar{\bar{\alpha}}\nabla\phi$：

$$ w\,\nabla\cdot(\bar{\bar{\alpha}}\nabla\phi) = \nabla\cdot\left(w\,\bar{\bar{\alpha}}\nabla\phi\right) - \nabla w\cdot(\bar{\bar{\alpha}}\nabla\phi) $$

因此：

$$ -\int_\Omega w\,\nabla\cdot(\bar{\bar{\alpha}}\nabla\phi)\,d\Omega = -\int_\Omega\nabla\cdot(w\,\bar{\bar{\alpha}}\nabla\phi)\,d\Omega + \int_\Omega \nabla w\cdot(\bar{\bar{\alpha}}\nabla\phi)\,d\Omega $$

2.3 对散度项应用 Gauss 散度定理

$$ \int_\Omega\nabla\cdot(w\,\bar{\bar{\alpha}}\nabla\phi)\,d\Omega = \oint_{\partial\Omega}w\,(\bar{\bar{\alpha}}\nabla\phi)\cdot\hat{n}\,d\Gamma $$

代入：

$$ -\int_\Omega w\,\nabla\cdot(\bar{\bar{\alpha}}\nabla\phi)\,d\Omega = -\oint_{\partial\Omega}w\,(\bar{\bar{\alpha}}\nabla\phi)\cdot\hat{n}\,d\Gamma + \int_\Omega \nabla w\cdot(\bar{\bar{\alpha}}\nabla\phi)\,d\Omega $$

2.4 展开梯度点积

$$ \nabla w\cdot(\bar{\bar{\alpha}}\nabla\phi) = \frac{\partial w}{\partial x}\alpha_x\frac{\partial \phi}{\partial x} + \frac{\partial w}{\partial y}\alpha_y\frac{\partial \phi}{\partial y} $$

2.5 代入边界条件

在 $\Gamma_2 = \partial\Omega$ 上，由 Robin 条件 (2)：

$$ (\bar{\bar{\alpha}}\nabla\phi)\cdot\hat{n} = q - \gamma\phi $$

代入边界积分：

$$ -\oint_{\partial\Omega}w\,(\bar{\bar{\alpha}}\nabla\phi)\cdot\hat{n}\,d\Gamma = -\oint_{\Gamma_2}w\,(q-\gamma\phi)\,d\Gamma $$$$ = -\oint_{\Gamma_2}w\,q\,d\Gamma + \oint_{\Gamma_2}w\,\gamma\,\phi\,d\Gamma $$

2.6 完整弱形式

将所有项合并到 (3) 中：

$$ \int_\Omega \nabla w\cdot(\bar{\bar{\alpha}}\nabla\phi)\,d\Omega + \int_\Omega w\,\beta\,\phi\,d\Omega + \oint_{\Gamma_2}w\,\gamma\,\phi\,d\Gamma = \int_\Omega w\,f\,d\Omega + \oint_{\Gamma_2}w\,q\,d\Gamma $$

展开写：

$$ \boxed{ \begin{aligned} &\int_\Omega\left(\alpha_x\frac{\partial w}{\partial x}\frac{\partial \phi}{\partial x}+\alpha_y\frac{\partial w}{\partial y}\frac{\partial \phi}{\partial y}\right)d\Omega \\ &\quad + \int_\Omega \beta\,w\,\phi\,d\Omega \\ &\quad + \oint_{\Gamma_2}\gamma\,w\,\phi\,d\Gamma \\ &= \int_\Omega f\,w\,d\Omega + \oint_{\Gamma_2}q\,w\,d\Gamma \end{aligned} } \tag{4} $$

各项的物理含义：

项表达式来源对应矩阵
刚度项$\int_\Omega \alpha_x\frac{\partial w}{\partial x}\frac{\partial \phi}{\partial x}+\alpha_y\frac{\partial w}{\partial y}\frac{\partial \phi}{\partial y}\,d\Omega$PDE中的扩散/传播项$[K]$
质量项$\int_\Omega \beta\,w\,\phi\,d\Omega$PDE中的反应项$[M]$
Robin边界项$\oint_{\Gamma_2}\gamma\,w\,\phi\,d\Gamma$边界条件代入$[R]$
体源项$\int_\Omega f\,w\,d\Omega$PDE右端$\{b_f\}$
边界源项$\oint_{\Gamma_2}q\,w\,d\Gamma$边界条件右端$\{b_q\}$

与三维电磁FEM的类比：

  • 刚度项 $\leftrightarrow$ 旋度-旋度项 $\int(\nabla\times\mathbf{T})\cdot(\nabla\times\mathbf{E})\,dV$
  • 质量项 $\leftrightarrow$ $-k_0^2\int\mathbf{T}\cdot\bar{\bar{\varepsilon}}_r\mathbf{E}\,dV$
  • Robin边界项 $\leftrightarrow$ ABC面积分 $-jk_0\oint(\hat{n}\times\mathbf{T})\cdot(\hat{n}\times\mathbf{E})\,dS$

三、全局编号下的离散化

3.1 网格与全局节点编号

将 $\Omega$ 剖分为 $M$ 个三角形单元，共有 $N$ 个全局节点，编号为 $I = 1, 2, \ldots, N$。

每个全局节点 $I$ 的坐标为 $(x_I, y_I)$。

3.2 全局基函数 $\mathcal{N}_I(x,y)$

每个全局节点 $I$ 对应一个全局基函数 $\mathcal{N}_I(x,y)$，它是一个分片线性的标量函数，满足：

$$ \mathcal{N}_I(x_J, y_J) = \delta_{IJ}, \quad I, J = 1, \ldots, N $$

即在自己的节点上取值1，在其他所有节点上取值0。

全局基函数的支撑域：

$$ \text{supp}(\mathcal{N}_I) = \bigcup_{\substack{e:\\\text{节点 }I\text{ 属于单元 }e}} \overline{\Omega^e} $$
           ╱╲
          ╱  ╲
         ╱ e₃ ╲
        ╱______╲
       ╱╲  I   ╱╲
      ╱  ╲    ╱  ╲
     ╱ e₂ ╲  ╱ e₄ ╲
    ╱______╲╱______╲
       ╱╲      ╱╲
      ╱  ╲    ╱  ╲
     ╱ e₁ ╲  ╱ e₅ ╲
    ╱______╲╱______╲

    N_I ≠ 0 在 e₁, e₂, ..., e₅ 内
    N_I = 0 在所有其他单元内
    N_I 在节点 I 处取值 1，像"帐篷"

在单元 $e$ 内：

$$ \mathcal{N}_I(x,y)\bigg|_{\Omega^e} = \begin{cases} N_k^{(e)}(x,y) & \text{若节点 } I \text{ 是单元 } e \text{ 的第 } k \text{ 个局部节点} \\[6pt] 0 & \text{若节点 } I \text{ 不是单元 } e \text{ 的任何节点} \end{cases} $$

与三维棱边元的对比：

  • 2D节点元：每个全局节点对应一个全局基函数，无方向符号
  • 3D棱边元：每个全局棱边对应一个全局基函数，有方向符号 $s_k^{(e)}$

节点元没有方向问题（标量），所以不需要符号修正，组装更简单。

3.3 场的全局展开

$$ \phi_h(x,y) = \sum_{I=1}^{N}\phi_I\,\mathcal{N}_I(x,y) $$

其中 $\phi_I = \phi_h(x_I, y_I)$ 是第 $I$ 个全局节点上的未知数（DOF）。

3.4 测试函数的选取（Galerkin法）

$$ w = \mathcal{N}_J, \quad J = 1, 2, \ldots, N $$

3.5 代入弱形式

将 $\phi_h = \sum_{I=1}^{N}\phi_I\,\mathcal{N}_I$ 和 $w = \mathcal{N}_J$ 代入弱形式 (4)：

左端第一项（刚度项）：

$$ \int_\Omega\left(\alpha_x\frac{\partial \mathcal{N}_J}{\partial x}\frac{\partial}{\partial x}\sum_I\phi_I\mathcal{N}_I + \alpha_y\frac{\partial \mathcal{N}_J}{\partial y}\frac{\partial}{\partial y}\sum_I\phi_I\mathcal{N}_I\right)d\Omega $$$$ = \sum_{I=1}^{N}\phi_I\int_\Omega\left(\alpha_x\frac{\partial \mathcal{N}_J}{\partial x}\frac{\partial \mathcal{N}_I}{\partial x} + \alpha_y\frac{\partial \mathcal{N}_J}{\partial y}\frac{\partial \mathcal{N}_I}{\partial y}\right)d\Omega $$$$ = \sum_{I=1}^{N} K_{JI}\,\phi_I $$

其中：

$$ K_{JI} = \int_\Omega\left(\alpha_x\frac{\partial \mathcal{N}_J}{\partial x}\frac{\partial \mathcal{N}_I}{\partial x} + \alpha_y\frac{\partial \mathcal{N}_J}{\partial y}\frac{\partial \mathcal{N}_I}{\partial y}\right)d\Omega $$

左端第二项（质量项）：

$$ \int_\Omega \beta\,\mathcal{N}_J\sum_I\phi_I\mathcal{N}_I\,d\Omega = \sum_{I=1}^{N}\phi_I\int_\Omega\beta\,\mathcal{N}_J\,\mathcal{N}_I\,d\Omega = \sum_{I=1}^{N}M_{JI}\,\phi_I $$

其中：

$$ M_{JI} = \int_\Omega\beta\,\mathcal{N}_J\,\mathcal{N}_I\,d\Omega $$

左端第三项（Robin边界项）：

$$ \oint_{\Gamma_2}\gamma\,\mathcal{N}_J\sum_I\phi_I\mathcal{N}_I\,d\Gamma = \sum_{I=1}^{N}\phi_I\oint_{\Gamma_2}\gamma\,\mathcal{N}_J\,\mathcal{N}_I\,d\Gamma = \sum_{I=1}^{N}R_{JI}\,\phi_I $$

其中：

$$ R_{JI} = \oint_{\Gamma_2}\gamma\,\mathcal{N}_J\,\mathcal{N}_I\,d\Gamma $$

右端项：

$$ b_J = \int_\Omega f\,\mathcal{N}_J\,d\Omega + \oint_{\Gamma_2}q\,\mathcal{N}_J\,d\Gamma $$

3.6 全局方程组

$$ \boxed{ \sum_{I=1}^{N}\left(K_{JI} + M_{JI} + R_{JI}\right)\phi_I = b_J, \quad J = 1, \ldots, N } $$

矩阵形式：

$$ \boxed{ \left([K] + [M] + [R]\right)\{\phi\} = \{b\} } $$$$ [A]\{\phi\} = \{b\}, \quad A_{JI} = K_{JI} + M_{JI} + R_{JI} $$

四、全域积分 → 单元积分之和

4.1 体积积分的分解

$$ K_{JI} = \int_\Omega\left(\alpha_x\frac{\partial \mathcal{N}_J}{\partial x}\frac{\partial \mathcal{N}_I}{\partial x} + \alpha_y\frac{\partial \mathcal{N}_J}{\partial y}\frac{\partial \mathcal{N}_I}{\partial y}\right)d\Omega = \sum_{e=1}^{M}\int_{\Omega^e}\left(\alpha_x\frac{\partial \mathcal{N}_J}{\partial x}\frac{\partial \mathcal{N}_I}{\partial x} + \alpha_y\frac{\partial \mathcal{N}_J}{\partial y}\frac{\partial \mathcal{N}_I}{\partial y}\right)d\Omega $$

4.2 支撑域筛选（与三维完全类比）

若节点 $I$ 不属于单元 $e$：

$$ \mathcal{N}_I(x,y)\bigg|_{\Omega^e} = 0, \quad \frac{\partial \mathcal{N}_I}{\partial x}\bigg|_{\Omega^e} = 0, \quad \frac{\partial \mathcal{N}_I}{\partial y}\bigg|_{\Omega^e} = 0 $$

被积函数中含 $\frac{\partial \mathcal{N}_I}{\partial x}$ 或 $\frac{\partial \mathcal{N}_I}{\partial y}$ 的因子为零，整个积分为零。

若节点 $J$ 不属于单元 $e$（但 $I$ 属于）：

同理，$\frac{\partial \mathcal{N}_J}{\partial x} = \frac{\partial \mathcal{N}_J}{\partial y} = 0$ 在 $\Omega^e$ 内，积分为零。

若节点 $I$ 和 $J$ 都属于单元 $e$：

设 $I$ 是单元 $e$ 的第 $i$ 个局部节点，$J$ 是第 $j$ 个局部节点。

$$ \mathcal{N}_I\bigg|_{\Omega^e} = N_i^{(e)}, \quad \mathcal{N}_J\bigg|_{\Omega^e} = N_j^{(e)} $$

积分非零。

4.3 结果

$$ \boxed{ \begin{aligned} K_{JI} &= \sum_{\substack{e=1\\I\in e,\;J\in e}}^{M} \underbrace{\int_{\Omega^e}\left(\alpha_x^{(e)}\frac{\partial N_j^{(e)}}{\partial x}\frac{\partial N_i^{(e)}}{\partial x} + \alpha_y^{(e)}\frac{\partial N_j^{(e)}}{\partial y}\frac{\partial N_i^{(e)}}{\partial y}\right)d\Omega}_{K^e_{ji}} \end{aligned} } $$

类似地：

$$ \boxed{ \begin{aligned} M_{JI} &= \sum_{\substack{e=1\\I\in e,\;J\in e}}^{M} \underbrace{\int_{\Omega^e}\beta^{(e)} N_j^{(e)} N_i^{(e)}\,d\Omega}_{M^e_{ji}} \end{aligned} } $$

对比三维棱边元：节点元没有符号 $s_i s_j$！

这是因为节点元的基函数是标量，没有方向问题。全局基函数在单元内的限制就是对应的局部基函数，不需要符号修正：

$$ \mathcal{N}_I\bigg|_{\Omega^e} = N_{i(I,e)}^{(e)} \quad \text{（无符号！）} $$

而三维棱边元需要：

$$ \boldsymbol{\mathcal{N}}_I\bigg|_{\Omega^e} = s_{i(I,e)}^{(e)}\,\mathbf{N}_{i(I,e)}^{(e)} \quad \text{（有符号！）} $$

4.4 边界积分的分解

边界 $\Gamma_2$ 由 $M_b$ 条边界线段组成：

$$ R_{JI} = \oint_{\Gamma_2}\gamma\,\mathcal{N}_J\,\mathcal{N}_I\,d\Gamma = \sum_{b=1}^{M_b}\int_{\Gamma^b}\gamma\,\mathcal{N}_J\,\mathcal{N}_I\,d\Gamma $$

每条边界线段 $\Gamma^b$ 有2个节点。只有当 $I$ 和 $J$ 都是 $\Gamma^b$ 的节点时，积分才非零。

$$ R_{JI} = \sum_{\substack{b=1\\I\in b,\;J\in b}}^{M_b}\underbrace{\int_{\Gamma^b}\gamma\,N_j^{(b)}\,N_i^{(b)}\,d\Gamma}_{R^b_{ji}} $$

类比三维：这里的边界线段积分 $\leftrightarrow$ 三维中的ABC三角形面积分

五、三角形单元的局部基函数

5.1 三角形单元定义

单元 $e$ 有3个节点，局部编号 $n_1, n_2, n_3$，坐标为 $(x_1^{(e)}, y_1^{(e)})$，$(x_2^{(e)}, y_2^{(e)})$，$(x_3^{(e)}, y_3^{(e)})$。

          n₃
          /\
         /  \
        /    \
       / Ω^e  \
      /________\
    n₁          n₂

5.2 面积坐标（重心坐标）

$$ \lambda_1 + \lambda_2 + \lambda_3 = 1, \quad \lambda_i \geq 0 $$$$ \begin{pmatrix}x\\y\\1\end{pmatrix} = \begin{pmatrix}x_1&x_2&x_3\\y_1&y_2&y_3\\1&1&1\end{pmatrix}\begin{pmatrix}\lambda_1\\\lambda_2\\\lambda_3\end{pmatrix} $$

反解：

$$ \lambda_i(x,y) = \frac{a_i + b_i x + c_i y}{2A^e} $$

其中 $A^e$ 是三角形面积：

$$ 2A^e = \det\begin{pmatrix}x_1&y_1&1\\x_2&y_2&1\\x_3&y_3&1\end{pmatrix} = (x_1-x_3)(y_2-y_3)-(x_2-x_3)(y_1-y_3) $$

系数：

$$ \begin{aligned} a_1 &= x_2 y_3 - x_3 y_2, \quad b_1 = y_2-y_3, \quad c_1 = x_3-x_2 \\ a_2 &= x_3 y_1 - x_1 y_3, \quad b_2 = y_3-y_1, \quad c_2 = x_1-x_3 \\ a_3 &= x_1 y_2 - x_2 y_1, \quad b_3 = y_1-y_2, \quad c_3 = x_2-x_1 \end{aligned} $$

5.3 线性基函数

三角形的线性基函数就是面积坐标本身：

$$ N_k^{(e)}(x,y) = \lambda_k(x,y), \quad k = 1, 2, 3 $$

关键性质：$N_k^{(e)}(x_j, y_j) = \delta_{kj}$

梯度（常数！）：

$$ \frac{\partial N_k^{(e)}}{\partial x} = \frac{b_k}{2A^e}, \quad \frac{\partial N_k^{(e)}}{\partial y} = \frac{c_k}{2A^e} $$

与三维四面体类似，线性基函数的梯度在单元内为常数。

5.4 面积坐标积分公式

$$ \int_{\Omega^e}\lambda_1^{p}\lambda_2^{q}\lambda_3^{r}\,d\Omega = \frac{p!\,q!\,r!}{(p+q+r+2)!}\cdot 2A^e $$

常用特殊情况：

$$ \int_{\Omega^e}\lambda_i\,d\Omega = \frac{A^e}{3} $$$$ \int_{\Omega^e}\lambda_i^2\,d\Omega = \frac{2!\cdot 0!\cdot 0!}{(2+2)!}\cdot 2A^e = \frac{A^e}{6} $$$$ \int_{\Omega^e}\lambda_i\lambda_j\,d\Omega\;(i\neq j) = \frac{1!\cdot 1!\cdot 0!}{(2+2)!}\cdot 2A^e = \frac{A^e}{12} $$

紧凑形式：

$$ \boxed{ \int_{\Omega^e}\lambda_i\lambda_j\,d\Omega = \frac{(1+\delta_{ij})}{12}\,A^e \cdot \frac{1}{1} = \begin{cases} A^e/6 & i=j\\ A^e/12 & i\neq j \end{cases} } $$

六、局部单元矩阵的计算

6.1 局部刚度矩阵 $[K^e]_{3\times 3}$

$$ K^e_{ji} = \int_{\Omega^e}\left(\alpha_x^{(e)}\frac{\partial N_j^{(e)}}{\partial x}\frac{\partial N_i^{(e)}}{\partial x} + \alpha_y^{(e)}\frac{\partial N_j^{(e)}}{\partial y}\frac{\partial N_i^{(e)}}{\partial y}\right)d\Omega $$

梯度为常数，$\alpha_x, \alpha_y$ 在单元内为常数，可以提出积分：

$$ K^e_{ji} = \left(\alpha_x^{(e)}\frac{b_j}{2A^e}\frac{b_i}{2A^e} + \alpha_y^{(e)}\frac{c_j}{2A^e}\frac{c_i}{2A^e}\right)\cdot A^e $$$$ \boxed{ K^e_{ji} = \frac{1}{4A^e}\left(\alpha_x^{(e)}\,b_j\,b_i + \alpha_y^{(e)}\,c_j\,c_i\right) } $$

完整的 $3\times 3$ 矩阵：

$$ [K^e] = \frac{1}{4A^e} \begin{pmatrix} \alpha_x b_1 b_1+\alpha_y c_1 c_1 & \alpha_x b_1 b_2+\alpha_y c_1 c_2 & \alpha_x b_1 b_3+\alpha_y c_1 c_3 \\ \alpha_x b_2 b_1+\alpha_y c_2 c_1 & \alpha_x b_2 b_2+\alpha_y c_2 c_2 & \alpha_x b_2 b_3+\alpha_y c_2 c_3 \\ \alpha_x b_3 b_1+\alpha_y c_3 c_1 & \alpha_x b_3 b_2+\alpha_y c_3 c_2 & \alpha_x b_3 b_3+\alpha_y c_3 c_3 \end{pmatrix} $$

6.1 局部刚度矩阵（续）

向量形式（更优雅）：

定义 $\mathbf{g}_k = (b_k, c_k)^T$，$\bar{\bar{D}} = \text{diag}(\alpha_x^{(e)}, \alpha_y^{(e)})$，则：

$$ \boxed{ K^e_{ji} = \frac{1}{4A^e}\,\mathbf{g}_j^T\,\bar{\bar{D}}\,\mathbf{g}_i = \frac{1}{4A^e}\begin{pmatrix}b_j & c_j\end{pmatrix}\begin{pmatrix}\alpha_x & 0\\0 & \alpha_y\end{pmatrix}\begin{pmatrix}b_i\\c_i\end{pmatrix} } $$

类比三维：三维棱边元的刚度矩阵 $S^e_{ji} = 4(\nabla\lambda_{a_j}\times\nabla\lambda_{b_j})\cdot(\nabla\lambda_{a_i}\times\nabla\lambda_{b_i})\cdot V^e$，也是由梯度的"内积"乘以体积。二维中 $\frac{1}{4A^e}\mathbf{g}_j^T\bar{\bar{D}}\mathbf{g}_i$ 扮演了完全相同的角色。

6.2 局部质量矩阵 $[M^e]_{3\times 3}$

$$ M^e_{ji} = \int_{\Omega^e}\beta^{(e)}\,N_j^{(e)}\,N_i^{(e)}\,d\Omega = \beta^{(e)}\int_{\Omega^e}\lambda_j\,\lambda_i\,d\Omega $$

代入面积坐标积分公式：

$$ \boxed{ \begin{aligned} M^e_{ji} &= \beta^{(e)}\cdot\frac{(1+\delta_{ji})}{12}\cdot A^e \cdot 2 \\ &= \beta^{(e)}\,A^e\,\frac{(1+\delta_{ji})}{12} \end{aligned} } $$

修正：重新代入面积坐标积分公式。

$\int_{\Omega^e}\lambda_j\lambda_i\,d\Omega = \frac{1!\cdot 1!\cdot 0!}{(1+1+0+2)!}\cdot 2A^e = \frac{1}{24}\cdot 2A^e = \frac{A^e}{12}$，当 $i\neq j$

$\int_{\Omega^e}\lambda_i^2\,d\Omega = \frac{2!\cdot 0!\cdot 0!}{(2+0+0+2)!}\cdot 2A^e = \frac{2}{24}\cdot 2A^e = \frac{A^e}{6}$，当 $i=j$

$$ M^e_{ji} = \begin{cases} \beta^{(e)}\,\frac{A^e}{6} & j = i \\[4pt] \beta^{(e)}\,\frac{A^e}{12} & j \neq i \end{cases} $$

完整的 $3\times 3$ 矩阵：

$$ \boxed{ [M^e] = \frac{\beta^{(e)}\,A^e}{12} \begin{pmatrix} 2 & 1 & 1 \\ 1 & 2 & 1 \\ 1 & 1 & 2 \end{pmatrix} } $$

类比三维：四面体中 $I_{pq} = (1+\delta_{pq})V^e/20$，与此完全对应，只是常数不同（三维是 $1/20$，二维是 $1/12$）。

6.3 局部右端矢量 $\{f^e\}_{3\times 1}$

$$ f^e_j = \int_{\Omega^e} f^{(e)}\,N_j^{(e)}\,d\Omega = f^{(e)}\int_{\Omega^e}\lambda_j\,d\Omega = f^{(e)}\cdot\frac{A^e}{3} $$

（假设 $f$ 在单元内为常数 $f^{(e)}$。）

$$ \boxed{ \{f^e\} = \frac{f^{(e)}\,A^e}{3}\begin{pmatrix}1\\1\\1\end{pmatrix} } $$

如果 $f$ 不是常数，用高斯积分或用节点值线性插值后解析积分：

$$ f^{(e)}(x,y) \approx \sum_{k=1}^{3}f_k\,\lambda_k $$$$ f^e_j = \int_{\Omega^e}\left(\sum_k f_k\lambda_k\right)\lambda_j\,d\Omega = \sum_k f_k\int_{\Omega^e}\lambda_k\lambda_j\,d\Omega $$$$ = \sum_k f_k\cdot\frac{(1+\delta_{kj})A^e}{12} = \frac{A^e}{12}(f_j + f_j + f_1 + f_2 + f_3 - f_j) = \frac{A^e}{12}(f_j + \sum_k f_k) $$

更清楚地写：

$$ f^e_j = \frac{A^e}{12}\left(2f_j + \sum_{k\neq j}f_k\right) = \frac{A^e}{12}\left(f_j + f_1 + f_2 + f_3\right) $$

七、边界单元矩阵（Robin条件）

7.1 边界线段定义

边界 $\Gamma_2$ 上的一条线段 $\Gamma^b$ 有2个节点，局部编号 $n_1^{(b)}, n_2^{(b)}$，长度为 $L^b$。

    n₁ ──────────── n₂
         Γ^b
         长度 L^b

7.2 一维线性基函数

在线段 $\Gamma^b$ 上，用参数 $t \in [0, 1]$：

$$ N_1^{(b)}(t) = 1 - t, \quad N_2^{(b)}(t) = t $$$$ d\Gamma = L^b\,dt $$

7.3 局部 Robin 矩阵 $[R^b]_{2\times 2}$

$$ R^b_{ji} = \int_{\Gamma^b}\gamma\,N_j^{(b)}\,N_i^{(b)}\,d\Gamma = \gamma^{(b)}\,L^b\int_0^1 N_j(t)\,N_i(t)\,dt $$

一维积分公式：

$$ \int_0^1 t^p(1-t)^q\,dt = \frac{p!\,q!}{(p+q+1)!} $$$$ \int_0^1 N_1^2\,dt = \int_0^1(1-t)^2\,dt = \frac{0!\cdot 2!}{3!} = \frac{1}{3} $$$$ \int_0^1 N_2^2\,dt = \int_0^1 t^2\,dt = \frac{2!\cdot 0!}{3!} = \frac{1}{3} $$$$ \int_0^1 N_1 N_2\,dt = \int_0^1 t(1-t)\,dt = \frac{1!\cdot 1!}{3!} = \frac{1}{6} $$$$ \boxed{ [R^b] = \frac{\gamma^{(b)}\,L^b}{6}\begin{pmatrix}2&1\\1&2\end{pmatrix} } $$

类比三维：三维 ABC 面矩阵 $[\tilde{B}^f]_{3\times 3}$ 用三角形面积坐标的 $(1+\delta_{mn})/12\cdot A^f$，这里是一维线积坐标的 $(1+\delta_{mn})/6\cdot L^b$。

维度递减：3D体→2D面→1D线，系数分母从 $20 \to 12 \to 6$，规律为 $\frac{(d+1)!}{2!} \cdot \frac{1}{1} = \frac{(d+1)!}{2}$，其中 $d$ 是积分域的维度。

7.4 局部边界源矢量 $\{q^b\}_{2\times 1}$

$$ q^b_j = \int_{\Gamma^b}q\,N_j^{(b)}\,d\Gamma = q^{(b)}\,L^b\int_0^1 N_j(t)\,dt = q^{(b)}\,L^b\cdot\frac{1}{2} $$$$ \boxed{ \{q^b\} = \frac{q^{(b)}\,L^b}{2}\begin{pmatrix}1\\1\end{pmatrix} } $$

如果 $q$ 在线段上线性变化，$q(t) = q_1(1-t) + q_2\,t$：

$$ q^b_j = L^b\int_0^1\left(q_1(1-t)+q_2\,t\right)N_j(t)\,dt $$$$ q^b_1 = L^b\left(q_1\int_0^1(1-t)^2\,dt + q_2\int_0^1 t(1-t)\,dt\right) = L^b\left(\frac{q_1}{3}+\frac{q_2}{6}\right) = \frac{L^b}{6}(2q_1+q_2) $$$$ q^b_2 = L^b\left(q_1\int_0^1 t(1-t)\,dt + q_2\int_0^1 t^2\,dt\right) = \frac{L^b}{6}(q_1+2q_2) $$$$ \{q^b\} = \frac{L^b}{6}\begin{pmatrix}2q_1+q_2\\q_1+2q_2\end{pmatrix} $$

八、最小算例：正方形上的 $2\times 2$ 网格

8.1 几何与网格

单位正方形 $\Omega = [0,1]\times[0,1]$，分为2个三角形。

    N3 (0,1)────────── N4 (1,1)
       |\              |
       |  \     Tri2   |
       |    \          |
       |      \        |
       |        \      |
       | Tri1     \    |
       |            \  |
       |              \|
    N1 (0,0)────────── N2 (1,0)

4个全局节点：

全局节点$x$$y$
$N_1$00
$N_2$10
$N_3$01
$N_4$11

2个三角形单元：

单元局部 $n_1$局部 $n_2$局部 $n_3$
Tri 1$N_1$$N_2$$N_3$
Tri 2$N_2$$N_4$$N_3$

4条边界线段（整个边界都是 $\Gamma_2$）：

边界段节点位置
$\Gamma^1$$N_1 \to N_2$下边
$\Gamma^2$$N_2 \to N_4$右边
$\Gamma^3$$N_4 \to N_3$上边
$\Gamma^4$$N_3 \to N_1$左边

8.2 全局-局部映射表

三角形单元的映射（节点元无符号！）：

         局部节点:    n₁    n₂    n₃
Tri 1 → 全局节点:    N1    N2    N3
Tri 2 → 全局节点:    N2    N4    N3

映射表 $\text{map}(e,k) = I$：

单元 $e$$k=1$$k=2$$k=3$
Tri 1$N_1$$N_2$$N_3$
Tri 2$N_2$$N_4$$N_3$

边界线段的映射：

边界段 $b$局部 $n_1^{(b)}$局部 $n_2^{(b)}$
$\Gamma^1$$N_1$$N_2$
$\Gamma^2$$N_2$$N_4$
$\Gamma^3$$N_4$$N_3$
$\Gamma^4$$N_3$$N_1$

与三维的关键区别：节点元没有方向符号 $s_k^{(e)}$。不需要 sign_map。组装时直接 $A_{JI} \mathrel{+}= K^e_{ji}$，不需要乘符号。

8.3 Tri 1 的局部矩阵计算

节点坐标：$(x_1,y_1)=(0,0)$，$(x_2,y_2)=(1,0)$，$(x_3,y_3)=(0,1)$

面积：

$$ 2A^{(1)} = (x_1-x_3)(y_2-y_3)-(x_2-x_3)(y_1-y_3) $$$$ = (0-0)(0-1)-(1-0)(0-1) = 0-(-1) = 1 $$$$ A^{(1)} = \frac{1}{2} $$

系数 $b_k, c_k$：

$$ b_1 = y_2-y_3 = 0-1 = -1, \quad c_1 = x_3-x_2 = 0-1 = -1 $$$$ b_2 = y_3-y_1 = 1-0 = 1, \quad c_2 = x_1-x_3 = 0-0 = 0 $$$$ b_3 = y_1-y_2 = 0-0 = 0, \quad c_3 = x_2-x_1 = 1-0 = 1 $$

正方形算例：$2\times 2$ 网格完整推导

9个节点 · 8个三角形 · 1个内部节点 · 8个边界节点

一、网格定义

1.1 几何与节点

单位正方形 $\Omega = [0,1]\times[0,1]$，用 $2\times 2$ 网格剖分。

9个全局节点：

    N7(0,1)────── N8(0.5,1)────── N9(1,1)
      |  \          |  \            |
      |    \  Tri6  |    \  Tri8   |
      |      \      |      \       |
      | Tri5   \    | Tri7   \     |
      |          \  |          \   |
    N4(0,0.5)──── N5(0.5,0.5)──── N6(1,0.5)
      |  \          |  \            |
      |    \  Tri2  |    \  Tri4   |
      |      \      |      \       |
      | Tri1   \    | Tri3   \     |
      |          \  |          \   |
    N1(0,0)────── N2(0.5,0)────── N3(1,0)
节点$x$$y$类型
$N_1$00边界（角点）
$N_2$0.50边界（底边）
$N_3$10边界（角点）
$N_4$00.5边界（左边）
$N_5$0.50.5内部
$N_6$10.5边界（右边）
$N_7$01边界（角点）
$N_8$0.51边界（上边）
$N_9$11边界（角点）

1.2 三角形单元

8个三角形，每个方格沿对角线分为2个：

单元$n_1$（局部）$n_2$（局部）$n_3$（局部）位置
Tri 1$N_1$$N_2$$N_5$左下↙
Tri 2$N_1$$N_5$$N_4$左下↗
Tri 3$N_2$$N_3$$N_6$右下↙
Tri 4$N_2$$N_6$$N_5$右下↗
Tri 5$N_4$$N_5$$N_8$左上↙
Tri 6$N_4$$N_8$$N_7$左上↗
Tri 7$N_5$$N_6$$N_9$右上↙
Tri 8$N_5$$N_9$$N_8$右上↗

1.3 边界线段

8条边界线段，覆盖正方形的4条边：

边界段 $b$局部 $n_1^{(b)}$局部 $n_2^{(b)}$位置
$\Gamma^1$$N_1$$N_2$底边左半
$\Gamma^2$$N_2$$N_3$底边右半
$\Gamma^3$$N_3$$N_6$右边下半
$\Gamma^4$$N_6$$N_9$右边上半
$\Gamma^5$$N_9$$N_8$上边右半
$\Gamma^6$$N_8$$N_7$上边左半
$\Gamma^7$$N_7$$N_4$左边上半
$\Gamma^8$$N_4$$N_1$左边下半

1.4 全局-局部映射表

体积单元映射 $\text{map}(e,k)$：

单元 $e$$k=1$$k=2$$k=3$
Tri 1125
Tri 2154
Tri 3236
Tri 4265
Tri 5458
Tri 6487
Tri 7569
Tri 8598

边界段映射 $\text{map}_b(b,k)$：

边界段 $b$$k=1$$k=2$
$\Gamma^1$12
$\Gamma^2$23
$\Gamma^3$36
$\Gamma^4$69
$\Gamma^5$98
$\Gamma^6$87
$\Gamma^7$74
$\Gamma^8$41

与三维的对比：节点元不需要 sign_map，因为标量基函数没有方向性。这是比三维棱边元简单的地方。

二、内部节点与边界节点——全局基函数的支撑域

2.1 内部节点 $N_5$ 的全局基函数

$N_5$ 属于所有8个三角形，因此：

$$ \text{supp}(\mathcal{N}_5) = \Omega^{(1)}\cup\Omega^{(2)}\cup\Omega^{(3)}\cup\Omega^{(4)}\cup\Omega^{(5)}\cup\Omega^{(6)}\cup\Omega^{(7)}\cup\Omega^{(8)} = \Omega $$
    N7 ─────── N8 ─────── N9
      |\        |\         |
      |  \      |  \       |
      |    \    |    \     |
      |  全部包含 N5 的 8 个  |
      |  三角形构成帐篷函数   |
      |        \|        \ |
    N4 ─────── N5 ─────── N6
      |\        |\         |
      |  \      |  \       |
      |    \    |    \     |
      |      \  |      \   |
      |        \|        \ |
    N1 ─────── N2 ─────── N3
    
    N_5 在中心取值1，在所有相邻节点取值0
    形成一个"帐篷"形状

$\mathcal{N}_5$ 在每个单元内的表达：

| 单元 | $\mathcal{N}_5\big|_{\Omega^e}$ | 局部编号 | |——|——————————-|———| | Tri 1 | $N_3^{(1)} = \lambda_3^{(1)}$ | 第3个局部节点 | | Tri 2 | $N_2^{(2)} = \lambda_2^{(2)}$ | 第2个局部节点 | | Tri 3 | 不属于 → $0$ | — | | Tri 4 | $N_3^{(4)} = \lambda_3^{(4)}$ | 第3个局部节点 | | Tri 5 | $N_2^{(5)} = \lambda_2^{(5)}$ | 第2个局部节点 | | Tri 6 | 不属于 → $0$ | — | | Tri 7 | $N_1^{(7)} = \lambda_1^{(7)}$ | 第1个局部节点 | | Tri 8 | $N_1^{(8)} = \lambda_1^{(8)}$ | 第1个局部节点 |

等等——让我重新检查。$N_5$ 在每个单元中的映射：

从映射表看，$N_5$ 出现在：

  • Tri 1: $k=3$（$\text{map}(1,3)=5$）→ $\mathcal{N}_5|_{\Omega^1} = N_3^{(1)}$
  • Tri 2: $k=2$（$\text{map}(2,2)=5$）→ $\mathcal{N}_5|_{\Omega^2} = N_2^{(2)}$
  • Tri 4: $k=3$（$\text{map}(4,3)=5$）→ $\mathcal{N}_5|_{\Omega^4} = N_3^{(4)}$
  • Tri 5: $k=2$（$\text{map}(5,2)=5$）→ $\mathcal{N}_5|_{\Omega^5} = N_2^{(5)}$
  • Tri 7: $k=1$（$\text{map}(7,1)=5$）→ $\mathcal{N}_5|_{\Omega^7} = N_1^{(7)}$
  • Tri 8: $k=1$（$\text{map}(8,1)=5$）→ $\mathcal{N}_5|_{\Omega^8} = N_1^{(8)}$

$N_5$ 不出现在 Tri 3（节点2,3,6）和 Tri 6（节点4,8,7）中。

所以 $\mathcal{N}_5$ 的支撑域是 6 个三角形，不是全部 8 个。

2.2 边界节点 $N_2$ 的全局基函数

$N_2$ 出现在：

  • Tri 1: $k=2$
  • Tri 3: $k=1$
  • Tri 4: $k=1$

（3个三角形）

$$ \text{supp}(\mathcal{N}_2) = \Omega^{(1)}\cup\Omega^{(3)}\cup\Omega^{(4)} $$

此外 $N_2$ 还在2条边界线段上：$\Gamma^1$（$k=2$）和 $\Gamma^2$（$k=1$）。

2.3 角点 $N_1$ 的全局基函数

$N_1$ 出现在：

  • Tri 1: $k=1$
  • Tri 2: $k=1$

（2个三角形）

$$ \text{supp}(\mathcal{N}_1) = \Omega^{(1)}\cup\Omega^{(2)} $$

边界线段：$\Gamma^1$（$k=1$）和 $\Gamma^8$（$k=2$）。

2.4 全局基函数在单元中的完整列表

$N_5$ 在各单元中的角色（内部节点的完整图景）：

$$ \mathcal{N}_5(x,y)\bigg|_{\Omega^e} = \begin{cases} \lambda_3^{(1)}(x,y) & e=1 \\ \lambda_2^{(2)}(x,y) & e=2 \\ 0 & e=3 \\ \lambda_3^{(4)}(x,y) & e=4 \\ \lambda_2^{(5)}(x,y) & e=5 \\ 0 & e=6 \\ \lambda_1^{(7)}(x,y) & e=7 \\ \lambda_1^{(8)}(x,y) & e=8 \end{cases} $$

关键观察：同一个全局基函数 $\mathcal{N}_5$ 在不同单元中对应不同的局部编号（有时是 $\lambda_1$，有时是 $\lambda_2$，有时是 $\lambda_3$）。这就是为什么需要映射表。

三、局部矩阵的完整数值计算

3.1 设定材料参数

取常数参数以便验证：

$$ \alpha_x = 1, \quad \alpha_y = 2, \quad \beta = 3+j, \quad \gamma = j, \quad f = 1, \quad q = 0 $$

$\beta$ 为复数体现 Helmholtz 方程的特性（如 $\beta = k^2$ 可以是复数）。$\gamma = j$ 对应类似 ABC 的阻抗边界条件。

3.2 Tri 1 的完整计算

节点：$n_1=N_1(0,0)$，$n_2=N_2(0.5,0)$，$n_3=N_5(0.5,0.5)$

面积：

$$ 2A^{(1)} = (x_1-x_3)(y_2-y_3)-(x_2-x_3)(y_1-y_3) $$$$ = (0-0.5)(0-0.5)-(0.5-0.5)(0-0.5) = (-0.5)(-0.5)-0 = 0.25 $$$$ A^{(1)} = 0.125 $$

系数 $b_k, c_k$：

$$ b_1 = y_2-y_3 = 0-0.5 = -0.5, \quad c_1 = x_3-x_2 = 0.5-0.5 = 0 $$$$ b_2 = y_3-y_1 = 0.5-0 = 0.5, \quad c_2 = x_1-x_3 = 0-0.5 = -0.5 $$$$ b_3 = y_1-y_2 = 0-0 = 0, \quad c_3 = x_2-x_1 = 0.5-0 = 0.5 $$

验证：$b_1+b_2+b_3 = -0.5+0.5+0 = 0$ ✓，$c_1+c_2+c_3 = 0-0.5+0.5 = 0$ ✓

刚度矩阵 $[K^{(1)}]$：

$$ K^{(1)}_{ji} = \frac{1}{4\times 0.125}\left(1\cdot b_j b_i + 2\cdot c_j c_i\right) = 2\left(b_j b_i + 2\,c_j c_i\right) $$

逐元素计算：

$$ K^{(1)}_{11} = 2[(-0.5)(-0.5)+2(0)(0)] = 2[0.25] = 0.5 $$$$ K^{(1)}_{12} = 2[(-0.5)(0.5)+2(0)(-0.5)] = 2[-0.25] = -0.5 $$$$ K^{(1)}_{13} = 2[(-0.5)(0)+2(0)(0.5)] = 2[0] = 0 $$$$ K^{(1)}_{22} = 2[(0.5)(0.5)+2(-0.5)(-0.5)] = 2[0.25+0.5] = 1.5 $$$$ K^{(1)}_{23} = 2[(0.5)(0)+2(-0.5)(0.5)] = 2[-0.5] = -1.0 $$$$ K^{(1)}_{33} = 2[(0)(0)+2(0.5)(0.5)] = 2[0.5] = 1.0 $$$$ [K^{(1)}] = \begin{pmatrix} 0.5 & -0.5 & 0 \\ -0.5 & 1.5 & -1.0 \\ 0 & -1.0 & 1.0 \end{pmatrix} $$

验证：对称 ✓；行和 = $(0, 0, 0)$ ✓（常数函数的梯度为零）

质量矩阵 $[M^{(1)}]$：

$$ [M^{(1)}] = \frac{(3+j)\times 0.125}{12}\begin{pmatrix}2&1&1\\1&2&1\\1&1&2\end{pmatrix} = \frac{3+j}{96}\begin{pmatrix}2&1&1\\1&2&1\\1&1&2\end{pmatrix} $$

数值（保留精确分数）：

$$ \frac{3+j}{96} = 0.03125 + 0.01042j $$$$ [M^{(1)}] = (0.03125+0.01042j)\begin{pmatrix}2&1&1\\1&2&1\\1&1&2\end{pmatrix} $$

右端矢量 $\{f^{(1)}\}$：

$$ \{f^{(1)}\} = \frac{1\times 0.125}{3}\begin{pmatrix}1\\1\\1\end{pmatrix} = \frac{1}{24}\begin{pmatrix}1\\1\\1\end{pmatrix} \approx \begin{pmatrix}0.04167\\0.04167\\0.04167\end{pmatrix} $$

3.3 Tri 2 的完整计算

节点：$n_1=N_1(0,0)$，$n_2=N_5(0.5,0.5)$，$n_3=N_4(0,0.5)$

面积：

$$ 2A^{(2)} = (0-0)(0.5-0.5)-(0.5-0)(0-0.5) = 0-(-0.25) = 0.25 $$$$ A^{(2)} = 0.125 $$

系数：

$$ b_1 = 0.5-0.5=0, \quad c_1 = 0-0.5=-0.5 $$$$ b_2 = 0.5-0=0.5, \quad c_2 = 0-0.5=-0.5 $$

等等——让我更仔细地计算：

$n_1=N_1(0,0)$，$n_2=N_5(0.5,0.5)$，$n_3=N_4(0,0.5)$

$$ b_1 = y_2-y_3 = 0.5-0.5 = 0 $$$$ c_1 = x_3-x_2 = 0-0.5 = -0.5 $$$$ b_2 = y_3-y_1 = 0.5-0 = 0.5 $$$$ c_2 = x_1-x_3 = 0-0 = 0 $$$$ b_3 = y_1-y_2 = 0-0.5 = -0.5 $$$$ c_3 = x_2-x_1 = 0.5-0 = 0.5 $$

验证：$b_1+b_2+b_3 = 0+0.5-0.5=0$ ✓，$c_1+c_2+c_3 = -0.5+0+0.5=0$ ✓

刚度矩阵：

$$ K^{(2)}_{ji} = 2(b_j b_i + 2\,c_j c_i) $$$$ K^{(2)}_{11} = 2[0+2(0.25)] = 1.0 $$$$ K^{(2)}_{12} = 2[0+2(-0.5)(0)] = 0 $$$$ K^{(2)}_{13} = 2[0+2(-0.5)(0.5)] = -1.0 $$$$ K^{(2)}_{22} = 2[0.25+0] = 0.5 $$$$ K^{(2)}_{23} = 2[(0.5)(-0.5)+0] = -0.5 $$$$ K^{(2)}_{33} = 2[0.25+2(0.25)] = 1.5 $$$$ [K^{(2)}] = \begin{pmatrix} 1.0 & 0 & -1.0 \\ 0 & 0.5 & -0.5 \\ -1.0 & -0.5 & 1.5 \end{pmatrix} $$

验证：对称 ✓；行和 = $(0,0,0)$ ✓

质量矩阵和右端项与 Tri 1 结构相同（$A^{(2)}=0.125$）：

$$ [M^{(2)}] = \frac{3+j}{96}\begin{pmatrix}2&1&1\\1&2&1\\1&1&2\end{pmatrix}, \quad \{f^{(2)}\} = \frac{1}{24}\begin{pmatrix}1\\1\\1\end{pmatrix} $$

3.4 对称性观察：所有8个三角形

由于网格的对称性，8个三角形分为两种类型：

类型 A（下三角，↙方向）：Tri 1, 3, 5, 7

这些三角形形状相同（直角在左下），$b_k, c_k$ 的值相同（可能差一个排列），所以刚度矩阵结构相同。

类型 B（上三角，↗方向）：Tri 2, 4, 6, 8

同理。

因此只需计算两种局部矩阵。

让我验证 Tri 3 和 Tri 7 与 Tri 1 结构相同：

Tri 3：$n_1=N_2(0.5,0)$，$n_2=N_3(1,0)$，$n_3=N_6(1,0.5)$

$$ b_1=0-0.5=-0.5,\; c_1=1-1=0 $$$$ b_2=0.5-0=0.5,\; c_2=0.5-1=-0.5 $$$$ b_3=0-0=0,\; c_3=1-0.5=0.5 $$

与 Tri 1 的 $(b_k,c_k)$ 完全相同！✓

Tri 4：$n_1=N_2(0.5,0)$，$n_2=N_6(1,0.5)$，$n_3=N_5(0.5,0.5)$

$$ b_1=0.5-0.5=0,\; c_1=0.5-1=-0.5 $$$$ b_2=0.5-0=0.5,\; c_2=0.5-0.5=0 $$$$ b_3=0-0.5=-0.5,\; c_3=1-0.5=0.5 $$

与 Tri 2 的 $(b_k,c_k)$ 完全相同！✓（只是局部节点排列不同，但 $K^e$ 的矩阵结构相同）

等等——让我仔细看。Tri 2 的 $(b_k,c_k)$：

Tri 2: $(0,-0.5)$, $(0.5,0)$, $(-0.5,0.5)$

Tri 4: $(0,-0.5)$, $(0.5,0)$, $(-0.5,0.5)$

完全相同！✓

结论：

$$ [K^A] \equiv [K^{(1)}] = [K^{(3)}] = [K^{(5)}] = [K^{(7)}] = \begin{pmatrix}0.5&-0.5&0\\-0.5&1.5&-1.0\\0&-1.0&1.0\end{pmatrix} $$$$ [K^B] \equiv [K^{(2)}] = [K^{(4)}] = [K^{(6)}] = [K^{(8)}] = \begin{pmatrix}1.0&0&-1.0\\0&0.5&-0.5\\-1.0&-0.5&1.5\end{pmatrix} $$

所有三角形面积相同 $A^e = 1/8$

四、边界线段矩阵的完整计算

4.1 所有边界段的参数

每条边界线段长度 $L^b = 0.5$（正方形边长1，每边分2段）。

Robin 矩阵（通用）：

$$ [R^b] = \frac{\gamma\,L^b}{6}\begin{pmatrix}2&1\\1&2\end{pmatrix} = \frac{j\times 0.5}{6}\begin{pmatrix}2&1\\1&2\end{pmatrix} = \frac{j}{12}\begin{pmatrix}2&1\\1&2\end{pmatrix} $$

边界源矢量（$q=0$）：

$$ \{q^b\} = \frac{0\times 0.5}{2}\begin{pmatrix}1\\1\end{pmatrix} = \begin{pmatrix}0\\0\end{pmatrix} $$

所有8条边界段的 $[R^b]$ 和 $\{q^b\}$ 结构完全相同，只是映射到不同的全局节点。

4.2 边界段映射汇总

边界段 $b$全局节点 $(I_1, I_2)$$[R^b]$ 贡献位置
$\Gamma^1$$(1, 2)$$A_{11}, A_{12}, A_{21}, A_{22}$
$\Gamma^2$$(2, 3)$$A_{22}, A_{23}, A_{32}, A_{33}$
$\Gamma^3$$(3, 6)$$A_{33}, A_{36}, A_{63}, A_{66}$
$\Gamma^4$$(6, 9)$$A_{66}, A_{69}, A_{96}, A_{99}$
$\Gamma^5$$(9, 8)$$A_{99}, A_{98}, A_{89}, A_{88}$
$\Gamma^6$$(8, 7)$$A_{88}, A_{87}, A_{78}, A_{77}$
$\Gamma^7$$(7, 4)$$A_{77}, A_{74}, A_{47}, A_{44}$
$\Gamma^8$$(4, 1)$$A_{44}, A_{41}, A_{14}, A_{11}$

关键观察：内部节点 $N_5$ 不出现在任何边界段中，所以 Robin 条件对第5行和第5列没有任何贡献。

五、全局矩阵的组装——完整过程

5.1 组装公式（节点元版本）

$$ A_{JI} = \underbrace{\sum_{\substack{e:\\I\in e,\,J\in e}} K^e_{j(J,e),\,i(I,e)}}_{\text{刚度}} + \underbrace{\sum_{\substack{e:\\I\in e,\,J\in e}} M^e_{j(J,e),\,i(I,e)}}_{\text{质量}} + \underbrace{\sum_{\substack{b:\\I\in b,\,J\in b}} R^b_{j(J,b),\,i(I,b)}}_{\text{Robin边界}} $$$$ b_J = \sum_{\substack{e:\\J\in e}} f^e_{j(J,e)} + \sum_{\substack{b:\\J\in b}} q^b_{j(J,b)} $$

对比三维棱边元：三维版本有 $s_i^{(e)} s_j^{(e)}$ 符号因子，这里没有。这是节点元与棱边元的根本区别。

5.2 逐单元散布

初始化：$[A]_{9\times 9} = [0]$，$\{b\}_{9\times 1} = \{0\}$

定义单元矩阵 $[K^e + M^e]$ 为：

对于类型 A 单元（Tri 1, 3, 5, 7）：

$$ [\hat{K}^A] \equiv [K^A] + [M^A] = \begin{pmatrix}0.5&-0.5&0\\-0.5&1.5&-1.0\\0&-1.0&1.0\end{pmatrix} + \frac{3+j}{96}\begin{pmatrix}2&1&1\\1&2&1\\1&1&2\end{pmatrix} $$

为简化组装演示，我先只用符号表示，最后再代入数值。

记 $K^e_{ji} + M^e_{ji} \equiv \hat{K}^e_{ji}$。

Tri 1 的散布：$\text{map}(1,\cdot) = (1, 2, 5)$

局部 (j,i):    (1,1)  (1,2)  (1,3)  (2,1)  (2,2)  (2,3)  (3,1)  (3,2)  (3,3)
                 ↓      ↓      ↓      ↓      ↓      ↓      ↓      ↓      ↓
全局 (J,I):    (1,1)  (1,2)  (1,5)  (2,1)  (2,2)  (2,5)  (5,1)  (5,2)  (5,5)
$$ A_{1,1} \mathrel{+}= \hat{K}^{(1)}_{11}, \quad A_{1,2} \mathrel{+}= \hat{K}^{(1)}_{12}, \quad A_{1,5} \mathrel{+}= \hat{K}^{(1)}_{13} $$$$ A_{2,1} \mathrel{+}= \hat{K}^{(1)}_{21}, \quad A_{2,2} \mathrel{+}= \hat{K}^{(1)}_{22}, \quad A_{2,5} \mathrel{+}= \hat{K}^{(1)}_{23} $$$$ A_{5,1} \mathrel{+}= \hat{K}^{(1)}_{31}, \quad A_{5,2} \mathrel{+}= \hat{K}^{(1)}_{32}, \quad A_{5,5} \mathrel{+}= \hat{K}^{(1)}_{33} $$

Tri 2 的散布：$\text{map}(2,\cdot) = (1, 5, 4)$

$$ A_{1,1} \mathrel{+}= \hat{K}^{(2)}_{11}, \quad A_{1,5} \mathrel{+}= \hat{K}^{(2)}_{12}, \quad A_{1,4} \mathrel{+}= \hat{K}^{(2)}_{13} $$$$ A_{5,1} \mathrel{+}= \hat{K}^{(2)}_{21}, \quad A_{5,5} \mathrel{+}= \hat{K}^{(2)}_{22}, \quad A_{5,4} \mathrel{+}= \hat{K}^{(2)}_{23} $$$$ A_{4,1} \mathrel{+}= \hat{K}^{(2)}_{31}, \quad A_{4,5} \mathrel{+}= \hat{K}^{(2)}_{32}, \quad A_{4,4} \mathrel{+}= \hat{K}^{(2)}_{33} $$

Tri 3 的散布：$\text{map}(3,\cdot) = (2, 3, 6)$

$$ A_{2,2} \mathrel{+}= \hat{K}^{(3)}_{11}, \quad A_{2,3} \mathrel{+}= \hat{K}^{(3)}_{12}, \quad A_{2,6} \mathrel{+}= \hat{K}^{(3)}_{13} $$$$ A_{3,2} \mathrel{+}= \hat{K}^{(3)}_{21}, \quad A_{3,3} \mathrel{+}= \hat{K}^{(3)}_{22}, \quad A_{3,6} \mathrel{+}= \hat{K}^{(3)}_{23} $$$$ A_{6,2} \mathrel{+}= \hat{K}^{(3)}_{31}, \quad A_{6,3} \mathrel{+}= \hat{K}^{(3)}_{32}, \quad A_{6,6} \mathrel{+}= \hat{K}^{(3)}_{33} $$

Tri 4 的散布：$\text{map}(4,\cdot) = (2, 6, 5)$

$$ A_{2,2} \mathrel{+}= \hat{K}^{(4)}_{11}, \quad A_{2,6} \mathrel{+}= \hat{K}^{(4)}_{12}, \quad A_{2,5} \mathrel{+}= \hat{K}^{(4)}_{13} $$$$ A_{6,2} \mathrel{+}= \hat{K}^{(4)}_{21}, \quad A_{6,6} \mathrel{+}= \hat{K}^{(4)}_{22}, \quad A_{6,5} \mathrel{+}= \hat{K}^{(4)}_{23} $$$$ A_{5,2} \mathrel{+}= \hat{K}^{(4)}_{31}, \quad A_{5,6} \mathrel{+}= \hat{K}^{(4)}_{32}, \quad A_{5,5} \mathrel{+}= \hat{K}^{(4)}_{33} $$

Tri 5 的散布：$\text{map}(5,\cdot) = (4, 5, 8)$

$$ A_{4,4} \mathrel{+}= \hat{K}^{(5)}_{11}, \quad A_{4,5} \mathrel{+}= \hat{K}^{(5)}_{12}, \quad A_{4,8} \mathrel{+}= \hat{K}^{(5)}_{13} $$$$ A_{5,4} \mathrel{+}= \hat{K}^{(5)}_{21}, \quad A_{5,5} \mathrel{+}= \hat{K}^{(5)}_{22}, \quad A_{5,8} \mathrel{+}= \hat{K}^{(5)}_{23} $$$$ A_{8,4} \mathrel{+}= \hat{K}^{(5)}_{31}, \quad A_{8,5} \mathrel{+}= \hat{K}^{(5)}_{32}, \quad A_{8,8} \mathrel{+}= \hat{K}^{(5)}_{33} $$

Tri 6 的散布：$\text{map}(6,\cdot) = (4, 8, 7)$

$$ A_{4,4} \mathrel{+}= \hat{K}^{(6)}_{11}, \quad A_{4,8} \mathrel{+}= \hat{K}^{(6)}_{12}, \quad A_{4,7} \mathrel{+}= \hat{K}^{(6)}_{13} $$$$ A_{8,4} \mathrel{+}= \hat{K}^{(6)}_{21}, \quad A_{8,8} \mathrel{+}= \hat{K}^{(6)}_{22}, \quad A_{8,7} \mathrel{+}= \hat{K}^{(6)}_{23} $$$$ A_{7,4} \mathrel{+}= \hat{K}^{(6)}_{31}, \quad A_{7,8} \mathrel{+}= \hat{K}^{(6)}_{32}, \quad A_{7,7} \mathrel{+}= \hat{K}^{(6)}_{33} $$

Tri 7 的散布：$\text{map}(7,\cdot) = (5, 6, 9)$

$$ A_{5,5} \mathrel{+}= \hat{K}^{(7)}_{11}, \quad A_{5,6} \mathrel{+}= \hat{K}^{(7)}_{12}, \quad A_{5,9} \mathrel{+}= \hat{K}^{(7)}_{13} $$$$ A_{6,5} \mathrel{+}= \hat{K}^{(7)}_{21}, \quad A_{6,6} \mathrel{+}= \hat{K}^{(7)}_{22}, \quad A_{6,9} \mathrel{+}= \hat{K}^{(7)}_{23} $$$$ A_{9,5} \mathrel{+}= \hat{K}^{(7)}_{31}, \quad A_{9,6} \mathrel{+}= \hat{K}^{(7)}_{32}, \quad A_{9,9} \mathrel{+}= \hat{K}^{(7)}_{33} $$

Tri 8 的散布：$\text{map}(8,\cdot) = (5, 9, 8)$

$$ A_{5,5} \mathrel{+}= \hat{K}^{(8)}_{11}, \quad A_{5,9} \mathrel{+}= \hat{K}^{(8)}_{12}, \quad A_{5,8} \mathrel{+}= \hat{K}^{(8)}_{13} $$$$ A_{9,5} \mathrel{+}= \hat{K}^{(8)}_{21}, \quad A_{9,9} \mathrel{+}= \hat{K}^{(8)}_{22}, \quad A_{9,8} \mathrel{+}= \hat{K}^{(8)}_{23} $$$$ A_{8,5} \mathrel{+}= \hat{K}^{(8)}_{31}, \quad A_{8,9} \mathrel{+}= \hat{K}^{(8)}_{32}, \quad A_{8,8} \mathrel{+}= \hat{K}^{(8)}_{33} $$

5.3 追踪每个全局矩阵元素的来源

现在我把所有单元的贡献按全局位置汇总。这就是"从局部回到全局"的过程。

对角元素 $A_{5,5}$（内部节点）——来源最多：

$$ A_{5,5} = \hat{K}^{(1)}_{33} + \hat{K}^{(2)}_{22} + \hat{K}^{(4)}_{33} + \hat{K}^{(5)}_{22} + \hat{K}^{(7)}_{11} + \hat{K}^{(8)}_{11} $$

代入数值（仅刚度部分）：

$$ K_{5,5} = 1.0 + 0.5 + 1.0 + 0.5 + 0.5 + 1.0 = 4.5 $$

代入数值（仅质量部分），所有 $M^e_{kk} = \beta A^e/6 = (3+j)/48$：

$$ M_{5,5} = 6\times\frac{3+j}{48} = \frac{3+j}{8} $$

（因为 $N_5$ 出现在6个三角形中，每个贡献对角质量项。）

$R_{5,5} = 0$（$N_5$ 是内部节点，不在任何边界段上。）

$$ A_{5,5} = 4.5 + \frac{3+j}{8} + 0 = 4.5 + 0.375 + 0.125j = 4.875 + 0.125j $$

对角元素 $A_{1,1}$（角点）——来源最少：

$$ A_{1,1} = \hat{K}^{(1)}_{11} + \hat{K}^{(2)}_{11} + R^{\Gamma^1}_{11} + R^{\Gamma^8}_{22} $$

其中 $R^{\Gamma^1}_{11} = \frac{j}{12}\cdot 2 = \frac{j}{6}$，$R^{\Gamma^8}_{22} = \frac{j}{12}\cdot 2 = \frac{j}{6}$

刚度：$K_{1,1} = 0.5 + 1.0 = 1.5$

质量：$M_{1,1} = 2\times\frac{3+j}{48} = \frac{3+j}{24}$

Robin：$R_{1,1} = \frac{j}{6} + \frac{j}{6} = \frac{j}{3}$

$$ A_{1,1} = 1.5 + \frac{3+j}{24} + \frac{j}{3} = 1.5 + 0.125 + \frac{j}{24} + \frac{j}{3} = 1.625 + \frac{9j}{24} = 1.625 + 0.375j $$

对角元素 $A_{2,2}$（边界中点）：

$N_2$ 出现在 Tri 1（$k=2$），Tri 3（$k=1$），Tri 4（$k=1$），共3个单元。

$$ K_{2,2} = K^{(1)}_{22} + K^{(3)}_{11} + K^{(4)}_{11} = 1.5 + 0.5 + 1.0 = 3.0 $$$$ M_{2,2} = 3\times\frac{3+j}{48} = \frac{3+j}{16} $$

$N_2$ 在边界段 $\Gamma^1$（$k=2$）和 $\Gamma^2$（$k=1$）上：

$$ R_{2,2} = R^{\Gamma^1}_{22} + R^{\Gamma^2}_{11} = \frac{j}{6} + \frac{j}{6} = \frac{j}{3} $$$$ A_{2,2} = 3.0 + \frac{3+j}{16} + \frac{j}{3} $$

非对角元素 $A_{5,2}$（内部-边界耦合）：

$N_5$ 和 $N_2$ 同时出现在哪些单元？

  • Tri 1：$N_2$ 是 $k=2$，$N_5$ 是 $k=3$ → 贡献 $\hat{K}^{(1)}_{32}$
  • Tri 4：$N_2$ 是 $k=1$，$N_5$ 是 $k=3$ → 贡献 $\hat{K}^{(4)}_{31}$
$$ K_{5,2} = K^{(1)}_{32} + K^{(4)}_{31} = (-1.0) + (-1.0) = -2.0 $$$$ M_{5,2} = M^{(1)}_{32} + M^{(4)}_{31} = \frac{3+j}{96} + \frac{3+j}{96} = \frac{3+j}{48} $$

$N_5$ 和 $N_2$ 不共享任何边界段，所以 $R_{5,2} = 0$。

$$ A_{5,2} = -2.0 + \frac{3+j}{48} $$

非对角元素 $A_{1,5}$：

$N_1$ 和 $N_5$ 同时出现在：

  • Tri 1：$N_1$ 是 $k=1$，$N_5$ 是 $k=3$ → 贡献 $\hat{K}^{(1)}_{13}$
  • Tri 2：$N_1$ 是 $k=1$，$N_5$ 是 $k=2$ → 贡献 $\hat{K}^{(2)}_{12}$
$$ K_{1,5} = K^{(1)}_{13} + K^{(2)}_{12} = 0 + 0 = 0 $$$$ M_{1,5} = M^{(1)}_{13} + M^{(2)}_{12} = \frac{3+j}{96} + \frac{3+j}{96} = \frac{3+j}{48} $$$$ A_{1,5} = \frac{3+j}{48} $$

有趣：刚度部分为零，只有质量项耦合 $N_1$ 和 $N_5$。

零元素 $A_{1,6}$：

$N_1$ 和 $N_6$ 不共享任何单元（$N_1$ 在 Tri 1,2；$N_6$ 在 Tri 3,4,7）。

由支撑域筛选：

$$ A_{1,6} = 0 $$

零元素 $A_{1,9}$：

$N_1$ 和 $N_9$ 相距最远，不共享任何单元。

$$ A_{1,9} = 0 $$

5.4 全局矩阵的稀疏结构

标记非零元素（$\bullet$）和零元素（$\cdot$）：

      N1   N2   N3   N4   N5   N6   N7   N8   N9
N1  [  •    •    ·    •    •    ·    ·    ·    ·  ]
N2  [  •    •    •    ·    •    •    ·    ·    ·  ]
N3  [  ·    •    •    ·    ·    •    ·    ·    ·  ]
N4  [  •    ·    ·    •    •    ·    •    •    ·  ]
N5  [  •    •    ·    •    •    •    ·    •    •  ]
N6  [  ·    •    •    ·    •    •    ·    ·    •  ]
N7  [  ·    ·    ·    •    ·    ·    •    •    ·  ]
N8  [  ·    ·    ·    •    •    ·    •    •    •  ]
N9  [  ·    ·    ·    ·    •    •    ·    •    •  ]

非零元素数：41个（$9\times 9 = 81$ 个中），稀疏率 ≈ 50%。

这里网格很小所以稀疏率不明显。实际大网格中（如 $100\times 100$ 网格，10201个节点），每行约5-7个非零元素，稀疏率 ≈ 0.05%。

5.5 完整的全局矩阵数值（仅刚度部分）

通过汇总所有8个单元的贡献：

$$ [K]_{global} = \begin{pmatrix} 1.5 & -0.5 & 0 & -1.0 & 0 & 0 & 0 & 0 & 0 \\ -0.5 & 3.0 & -0.5 & 0 & -2.0 & 0 & 0 & 0 & 0 \\ 0 & -0.5 & 1.5 & 0 & 0 & -1.0 & 0 & 0 & 0 \\ -1.0 & 0 & 0 & 3.0 & -2.0 & 0 & -0.5 & 0 & 0 \\ -0.5 0 & -2.0 & 0 & -2.0 & 8.0 & -2.0 & 0 & -2.0 & 0 \\ 待修正 0 & 0 & -1.0 & 0 & -2.0 & 3.0 & 0 & 0 & -0.5 \\ 待修正 0 & 0 & 0 & -0.5 & 0 & 0 & 1.5 & -0.5 & 0 \\ 待修正 0 & 0 & 0 & 0 & -2.0 & 0 & -0.5 & 3.0 & -0.5 \\ 待修正 0 & 0 & 0 & 0 & 0 & -0.5 & 0 & -0.5 & 1.5 \end{pmatrix} $$

让我仔细重新计算每个元素。

系统方法：列出每个全局节点对 $(J,I)$ 的所有单元贡献。

我先构造一个完整的贡献追踪表：

对于每个全局位置 $(J,I)$，列出所有贡献的单元和对应的局部位置。

$A_{4,4}$：$N_4$ 出现在 Tri 2($k=3$), Tri 5($k=1$), Tri 6($k=1$)

$$ K_{4,4} = K^{B}_{33} + K^{A}_{11} + K^{B}_{11} = 1.5 + 0.5 + 1.0 = 3.0 $$

$A_{4,5}$：$N_4$ 和 $N_5$ 共享 Tri 2（$N_4$:$k=3$, $N_5$:$k=2$）和 Tri 5（$N_4$:$k=1$, $N_5$:$k=2$）

$$ K_{4,5} = K^{B}_{32} + K^{A}_{12} = -0.5 + (-0.5) = -1.0 $$

$A_{4,1}$：$N_4$ 和 $N_1$ 共享 Tri 2（$N_4$:$k=3$, $N_1$:$k=1$）

$$ K_{4,1} = K^{B}_{31} = -1.0 $$

$A_{4,7}$：$N_4$ 和 $N_7$ 共享 Tri 6（$N_4$:$k=1$, $N_7$:$k=3$）

$$ K_{4,7} = K^{B}_{13} = -1.0 $$

$A_{4,8}$：$N_4$ 和 $N_8$ 共享 Tri 5（$N_4$:$k=1$, $N_8$:$k=3$）和 Tri 6（$N_4$:$k=1$, $N_8$:$k=2$）

$$ K_{4,8} = K^{A}_{13} + K^{B}_{12} = 0 + 0 = 0 $$

$A_{5,5}$（重新仔细计算）：

$N_5$ 出现在6个单元，列出每个的局部编号和对角刚度贡献：

单元类型$N_5$ 的局部编号 $k$$K^e_{kk}$
Tri 1A3$K^A_{33} = 1.0$
Tri 2B2$K^B_{22} = 0.5$
Tri 4B3$K^B_{33} = 1.5$
Tri 5A2$K^A_{22} = 1.5$
Tri 7A1$K^A_{11} = 0.5$
Tri 8B1$K^B_{11} = 1.0$
$$ K_{5,5} = 1.0 + 0.5 + 1.5 + 1.5 + 0.5 + 1.0 = 6.0 $$

$A_{5,2}$：$N_5$ 和 $N_2$ 共享 Tri 1 和 Tri 4

单元$N_5$ 局部$N_2$ 局部贡献
Tri 1 (A)32$K^A_{32}=-1.0$
Tri 4 (B)31$K^B_{31}=-1.0$
$$ K_{5,2} = -1.0 + (-1.0) = -2.0 $$

$A_{5,4}$：$N_5$ 和 $N_4$ 共享 Tri 2 和 Tri 5

单元$N_5$ 局部$N_4$ 局部贡献
Tri 2 (B)23$K^B_{23}=-0.5$
Tri 5 (A)21$K^A_{21}=-0.5$

五（续）：完整的全局矩阵组装

5.5 全局刚度矩阵——逐元素完整追踪

继续系统地计算每个全局矩阵元素。

回顾两种单元矩阵：

$$ [K^A] = \begin{pmatrix}0.5&-0.5&0\\-0.5&1.5&-1.0\\0&-1.0&1.0\end{pmatrix}, \quad [K^B] = \begin{pmatrix}1.0&0&-1.0\\0&0.5&-0.5\\-1.0&-0.5&1.5\end{pmatrix} $$

回顾映射表：

单元类型$k=1$$k=2$$k=3$
Tri 1A$N_1$$N_2$$N_5$
Tri 2B$N_1$$N_5$$N_4$
Tri 3A$N_2$$N_3$$N_6$
Tri 4B$N_2$$N_6$$N_5$
Tri 5A$N_4$$N_5$$N_8$
Tri 6B$N_4$$N_8$$N_7$
Tri 7A$N_5$$N_6$$N_9$
Tri 8B$N_5$$N_9$$N_8$

为每个全局节点，列出它在哪些单元中、局部编号是什么：

全局节点所在单元(局部编号)
$N_1$Tri1($k$=1), Tri2($k$=1)
$N_2$Tri1($k$=2), Tri3($k$=1), Tri4($k$=1)
$N_3$Tri3($k$=2)
$N_4$Tri2($k$=3), Tri5($k$=1), Tri6($k$=1)
$N_5$Tri1($k$=3), Tri2($k$=2), Tri4($k$=3), Tri5($k$=2), Tri7($k$=1), Tri8($k$=1)
$N_6$Tri3($k$=3), Tri4($k$=2), Tri7($k$=2)
$N_7$Tri6($k$=3)
$N_8$Tri5($k$=3), Tri6($k$=2), Tri8($k$=3)
$N_9$Tri7($k$=3), Tri8($k$=2)

判断两个节点是否共享单元的方法：取两个节点的"所在单元"集合的交集。

第1行：$J = N_1$

$N_1$ 在 Tri1($k$=1), Tri2($k$=1)

$K_{1,1}$：Tri1 ∩ Tri1, Tri2 ∩ Tri2

$$ K_{1,1} = K^A_{11} + K^B_{11} = 0.5 + 1.0 = 1.5 $$

$K_{1,2}$：$N_1$ 和 $N_2$ 共享 Tri1

单元$N_1$$N_2$贡献
Tri1(A)$k$=1$k$=2$K^A_{12} = -0.5$
$$ K_{1,2} = -0.5 $$

$K_{1,3}$：$N_1$(Tri1,Tri2) 与 $N_3$(Tri3) → 无共享单元

$$ K_{1,3} = 0 $$

$K_{1,4}$：$N_1$ 和 $N_4$ 共享 Tri2

单元$N_1$$N_4$贡献
Tri2(B)$k$=1$k$=3$K^B_{13} = -1.0$
$$ K_{1,4} = -1.0 $$

$K_{1,5}$：$N_1$ 和 $N_5$ 共享 Tri1, Tri2

单元$N_1$$N_5$贡献
Tri1(A)$k$=1$k$=3$K^A_{13} = 0$
Tri2(B)$k$=1$k$=2$K^B_{12} = 0$
$$ K_{1,5} = 0 + 0 = 0 $$

$K_{1,6}$：$N_1$(Tri1,2) 与 $N_6$(Tri3,4,7) → 无共享

$$ K_{1,6} = 0 $$

$K_{1,7} = K_{1,8} = K_{1,9} = 0$（均无共享单元）

第1行汇总：

$$ K_{1,\cdot} = (1.5,\;-0.5,\;0,\;-1.0,\;0,\;0,\;0,\;0,\;0) $$

验证：行和 $= 1.5-0.5-1.0 = 0$ ✓

第2行：$J = N_2$

$N_2$ 在 Tri1($k$=2), Tri3($k$=1), Tri4($k$=1)

$K_{2,1}$：共享 Tri1

$$ K_{2,1} = K^A_{21} = -0.5 $$

$K_{2,2}$：

$$ K_{2,2} = K^A_{22} + K^A_{11} + K^B_{11} = 1.5 + 0.5 + 1.0 = 3.0 $$

$K_{2,3}$：$N_2$ 和 $N_3$ 共享 Tri3

单元$N_2$$N_3$贡献
Tri3(A)$k$=1$k$=2$K^A_{12} = -0.5$
$$ K_{2,3} = -0.5 $$

$K_{2,4}$：$N_2$(Tri1,3,4) 与 $N_4$(Tri2,5,6) → 无共享

$$ K_{2,4} = 0 $$

$K_{2,5}$：共享 Tri1, Tri4

单元$N_2$$N_5$贡献
Tri1(A)$k$=2$k$=3$K^A_{23} = -1.0$
Tri4(B)$k$=1$k$=3$K^B_{13} = -1.0$
$$ K_{2,5} = -1.0 + (-1.0) = -2.0 $$

$K_{2,6}$：共享 Tri3, Tri4

单元$N_2$$N_6$贡献
Tri3(A)$k$=1$k$=3$K^A_{13} = 0$
Tri4(B)$k$=1$k$=2$K^B_{12} = 0$
$$ K_{2,6} = 0 + 0 = 0 $$

$K_{2,7} = K_{2,8} = K_{2,9} = 0$

第2行汇总：

$$ K_{2,\cdot} = (-0.5,\;3.0,\;-0.5,\;0,\;-2.0,\;0,\;0,\;0,\;0) $$

验证：行和 $= -0.5+3.0-0.5-2.0 = 0$ ✓

第3行：$J = N_3$

$N_3$ 仅在 Tri3($k$=2)

$K_{3,1} = 0$，$K_{3,4} = K_{3,5} = K_{3,7} = K_{3,8} = K_{3,9} = 0$

$K_{3,2}$：共享 Tri3

$$ K_{3,2} = K^A_{21} = -0.5 $$

$K_{3,3}$：

$$ K_{3,3} = K^A_{22} = 1.5 $$

$K_{3,6}$：共享 Tri3

$$ K_{3,6} = K^A_{23} = -1.0 $$

第3行汇总：

$$ K_{3,\cdot} = (0,\;-0.5,\;1.5,\;0,\;0,\;-1.0,\;0,\;0,\;0) $$

验证：行和 $= -0.5+1.5-1.0 = 0$ ✓

第4行：$J = N_4$

$N_4$ 在 Tri2($k$=3), Tri5($k$=1), Tri6($k$=1)

$K_{4,1}$：共享 Tri2

$$ K_{4,1} = K^B_{31} = -1.0 $$

$K_{4,2} = K_{4,3} = 0$

$K_{4,4}$：

$$ K_{4,4} = K^B_{33} + K^A_{11} + K^B_{11} = 1.5 + 0.5 + 1.0 = 3.0 $$

$K_{4,5}$：共享 Tri2, Tri5

单元$N_4$$N_5$贡献
Tri2(B)$k$=3$k$=2$K^B_{32} = -0.5$
Tri5(A)$k$=1$k$=2$K^A_{12} = -0.5$
$$ K_{4,5} = -0.5 + (-0.5) = -1.0 $$

$K_{4,6} = 0$

$K_{4,7}$：共享 Tri6

$$ K_{4,7} = K^B_{13} = -1.0 $$

$K_{4,8}$：共享 Tri5, Tri6

单元$N_4$$N_8$贡献
Tri5(A)$k$=1$k$=3$K^A_{13} = 0$
Tri6(B)$k$=1$k$=2$K^B_{12} = 0$
$$ K_{4,8} = 0 + 0 = 0 $$

$K_{4,9} = 0$

第4行汇总：

$$ K_{4,\cdot} = (-1.0,\;0,\;0,\;3.0,\;-1.0,\;0,\;-1.0,\;0,\;0) $$

验证：行和 $= -1.0+3.0-1.0-1.0 = 0$ ✓

第5行：$J = N_5$（内部节点——最复杂）

$N_5$ 在 Tri1($k$=3), Tri2($k$=2), Tri4($k$=3), Tri5($k$=2), Tri7($k$=1), Tri8($k$=1)

$K_{5,1}$：共享 Tri1, Tri2

单元$N_5$$N_1$贡献
Tri1(A)$k$=3$k$=1$K^A_{31} = 0$
Tri2(B)$k$=2$k$=1$K^B_{21} = 0$
$$ K_{5,1} = 0 $$

$K_{5,2}$：共享 Tri1, Tri4

单元$N_5$$N_2$贡献
Tri1(A)$k$=3$k$=2$K^A_{32} = -1.0$
Tri4(B)$k$=3$k$=1$K^B_{31} = -1.0$
$$ K_{5,2} = -2.0 $$

$K_{5,3} = 0$（无共享单元）

$K_{5,4}$：共享 Tri2, Tri5

单元$N_5$$N_4$贡献
Tri2(B)$k$=2$k$=3$K^B_{23} = -0.5$
Tri5(A)$k$=2$k$=1$K^A_{21} = -0.5$
$$ K_{5,4} = -1.0 $$

$K_{5,5}$：

单元$N_5$ 局部$K^e_{kk}$
Tri1(A)$k$=3$K^A_{33}=1.0$
Tri2(B)$k$=2$K^B_{22}=0.5$
Tri4(B)$k$=3$K^B_{33}=1.5$
Tri5(A)$k$=2$K^A_{22}=1.5$
Tri7(A)$k$=1$K^A_{11}=0.5$
Tri8(B)$k$=1$K^B_{11}=1.0$
$$ K_{5,5} = 1.0+0.5+1.5+1.5+0.5+1.0 = 6.0 $$

$K_{5,6}$：共享 Tri4, Tri7

单元$N_5$$N_6$贡献
Tri4(B)$k$=3$k$=2$K^B_{32} = -0.5$
Tri7(A)$k$=1$k$=2$K^A_{12} = -0.5$
$$ K_{5,6} = -1.0 $$

$K_{5,7} = 0$（$N_5$ 和 $N_7$ 无共享单元）

$K_{5,8}$：共享 Tri5, Tri8

单元$N_5$$N_8$贡献
Tri5(A)$k$=2$k$=3$K^A_{23} = -1.0$
Tri8(B)$k$=1$k$=3$K^B_{13} = -1.0$
$$ K_{5,8} = -2.0 $$

$K_{5,9}$：共享 Tri7, Tri8

单元$N_5$$N_9$贡献
Tri7(A)$k$=1$k$=3$K^A_{13} = 0$
Tri8(B)$k$=1$k$=2$K^B_{12} = 0$
$$ K_{5,9} = 0 $$

第5行汇总：

$$ K_{5,\cdot} = (0,\;-2.0,\;0,\;-1.0,\;6.0,\;-1.0,\;0,\;-2.0,\;0) $$

验证：行和 $= -2.0-1.0+6.0-1.0-2.0 = 0$ ✓

第6-9行（利用对称性快速给出）

由对称性 $K_{JI} = K_{IJ}$，可直接从前面结果读出列值。但为完整性，我逐行列出：

第6行：$N_6$ 在 Tri3($k$=3), Tri4($k$=2), Tri7($k$=2)

$$ K_{6,\cdot} = (0,\;0,\;-1.0,\;0,\;-1.0,\;3.0,\;0,\;0,\;-1.0) $$

验证：$-1.0-1.0+3.0-1.0=0$ ✓

让我仔细验证 $K_{6,9}$：$N_6$ 和 $N_9$ 共享 Tri7

单元$N_6$$N_9$贡献
Tri7(A)$k$=2$k$=3$K^A_{23} = -1.0$
$$ K_{6,9} = -1.0 \quad \checkmark $$

第7行：$N_7$ 仅在 Tri6($k$=3)

$$ K_{7,\cdot} = (0,\;0,\;0,\;-1.0,\;0,\;0,\;1.5,\;-0.5,\;0) $$

验证：$-1.0+1.5-0.5=0$ ✓

$K_{7,4}$：共享 Tri6 → $K^B_{31} = -1.0$ ✓

$K_{7,7}$：$K^B_{33} = 1.5$ ✓

$K_{7,8}$：共享 Tri6 → $K^B_{32} = -0.5$ ✓

第8行：$N_8$ 在 Tri5($k$=3), Tri6($k$=2), Tri8($k$=3)

$K_{8,4}$：共享 Tri5, Tri6

单元$N_8$$N_4$贡献
Tri5(A)$k$=3$k$=1$K^A_{31} = 0$
Tri6(B)$k$=2$k$=1$K^B_{21} = 0$
$$ K_{8,4} = 0 $$

$K_{8,5}$：共享 Tri5, Tri8

单元$N_8$$N_5$贡献
Tri5(A)$k$=3$k$=2$K^A_{32} = -1.0$
Tri8(B)$k$=3$k$=1$K^B_{31} = -1.0$
$$ K_{8,5} = -2.0 $$

$K_{8,7}$：共享 Tri6 → $K^B_{23} = -0.5$

$K_{8,8}$：

单元$N_8$ 局部$K^e_{kk}$
Tri5(A)$k$=3$K^A_{33}=1.0$
Tri6(B)$k$=2$K^B_{22}=0.5$
Tri8(B)$k$=3$K^B_{33}=1.5$
$$ K_{8,8} = 1.0+0.5+1.5 = 3.0 $$

$K_{8,9}$：共享 Tri8

单元$N_8$$N_9$贡献
Tri8(B)$k$=3$k$=2$K^B_{32} = -0.5$
$$ K_{8,9} = -0.5 $$

第8行汇总：

$$ K_{8,\cdot} = (0,\;0,\;0,\;0,\;-2.0,\;0,\;-0.5,\;3.0,\;-0.5) $$

验证：$-2.0-0.5+3.0-0.5=0$ ✓

第9行：$N_9$ 在 Tri7($k$=3), Tri8($k$=2)

$K_{9,5}$：共享 Tri7, Tri8

单元$N_9$$N_5$贡献
Tri7(A)$k$=3$k$=1$K^A_{31} = 0$
Tri8(B)$k$=2$k$=1$K^B_{21} = 0$
$$ K_{9,5} = 0 $$

$K_{9,6}$：共享 Tri7 → $K^A_{32} = -1.0$

$K_{9,8}$：共享 Tri8 → $K^B_{23} = -0.5$

$K_{9,9}$：

$$ K_{9,9} = K^A_{33} + K^B_{22} = 1.0 + 0.5 = 1.5 $$

第9行汇总：

$$ K_{9,\cdot} = (0,\;0,\;0,\;0,\;0,\;-1.0,\;0,\;-0.5,\;1.5) $$

验证：$-1.0-0.5+1.5=0$ ✓

5.6 完整的全局刚度矩阵

$$ \boxed{ [K]_{9\times 9} = \begin{pmatrix} 1.5 & -0.5 & 0 & -1.0 & 0 & 0 & 0 & 0 & 0 \\ -0.5 & 3.0 & -0.5 & 0 & -2.0 & 0 & 0 & 0 & 0 \\ 0 & -0.5 & 1.5 & 0 & 0 & -1.0 & 0 & 0 & 0 \\ -1.0 & 0 & 0 & 3.0 & -1.0 & 0 & -1.0 & 0 & 0 \\ 0 & -2.0 & 0 & -1.0 & 6.0 & -1.0 & 0 & -2.0 & 0 \\ 0 & 0 & -1.0 & 0 & -1.0 & 3.0 & 0 & 0 & -1.0 \\ 0 & 0 & 0 & -1.0 & 0 & 0 & 1.5 & -0.5 & 0 \\ 0 & 0 & 0 & 0 & -2.0 & 0 & -0.5 & 3.0 & -0.5 \\ 0 & 0 & 0 & 0 & 0 & -1.0 & 0 & -0.5 & 1.5 \end{pmatrix} } $$

全局验证：

  1. 对称性：$K_{JI} = K_{IJ}$ ✓（逐个检查）
  2. 每行和为零：所有9行 ✓（对应常数函数的梯度为零）
  3. 内部节点 $N_5$ 对角元最大（$6.0$），角点最小（$1.5$）✓

5.7 完整的全局质量矩阵

$$ [M]_{9\times 9} = \frac{\beta\,A^e}{12}\times\text{（组装后的矩阵）} $$

设 $\mu \equiv \frac{\beta\,A^e}{12} = \frac{(3+j)\times 0.125}{12} = \frac{3+j}{96}$

每个节点对 $(J,I)$ 的质量矩阵贡献为 $2\mu$（对角）或 $\mu$（非对角）乘以共享单元数。

内部节点 $N_5$ 的对角元：6个单元共享

$$ M_{5,5} = 6\times 2\mu = 12\mu = \frac{12(3+j)}{96} = \frac{3+j}{8} $$

边中点（如 $N_2$）的对角元：3个单元共享

$$ M_{2,2} = 3\times 2\mu = 6\mu = \frac{6(3+j)}{96} = \frac{3+j}{16} $$

角点（如 $N_1$）的对角元：2个单元共享

$$ M_{1,1} = 2\times 2\mu = 4\mu = \frac{4(3+j)}{96} = \frac{3+j}{24} $$

非对角元（共享单元数 $\times\mu$）：

$$ M_{5,2} = 2\mu, \quad M_{5,4} = 2\mu, \quad M_{5,6} = 2\mu, \quad M_{5,8} = 2\mu $$$$ M_{5,1} = 2\mu, \quad M_{5,9} = 2\mu $$$$ M_{1,2} = 1\mu, \quad M_{1,4} = 1\mu, \quad M_{1,5} = 2\mu $$

完整的全局质量矩阵：

$$ [M]_{9\times 9} = \mu\begin{pmatrix} 4 & 1 & 0 & 1 & 2 & 0 & 0 & 0 & 0 \\ 1 & 6 & 1 & 0 & 2 & 0 & 0 & 0 & 0 \\ 0 & 1 & 4 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 6 & 2 & 0 & 1 & 0 & 0 \\ 2 & 2 & 0 & 2 & 12 & 2 & 0 & 2 & 2 \\ 0 & 0 & 1 & 0 & 2 & 6 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & 4 & 1 & 0 \\ 0 & 0 & 0 & 0 & 2 & 0 & 1 & 6 & 1 \\ 0 & 0 & 0 & 0 & 2 & 1 & 0 & 1 & 4 \end{pmatrix} $$

等等，让我仔细检查 $M_{5,1}$ 和 $M_{5,9}$。

$M_{5,1}$：$N_5$ 和 $N_1$ 共享 Tri1, Tri2 → 2个共享单元，非对角 → $2\times\mu = 2\mu$ ✓

$M_{5,9}$：$N_5$ 和 $N_9$ 共享 Tri7, Tri8 → 2个共享单元 → $2\mu$ ✓

$M_{2,6}$：$N_2$ 和 $N_6$ 共享 Tri3, Tri4 → $2\mu$

我之前写成0，需要修正！

五（续）：质量矩阵的完整检查与修正

5.8 系统方法——列出每对节点的共享单元数

回顾每个全局节点所在的单元：

全局节点所在单元
$N_1$Tri1, Tri2
$N_2$Tri1, Tri3, Tri4
$N_3$Tri3
$N_4$Tri2, Tri5, Tri6
$N_5$Tri1, Tri2, Tri4, Tri5, Tri7, Tri8
$N_6$Tri3, Tri4, Tri7
$N_7$Tri6
$N_8$Tri5, Tri6, Tri8
$N_9$Tri7, Tri8

逐对计算共享单元数（取两个集合的交集大小）：

$(J,I)$$N_1$$N_2$$N_3$$N_4$$N_5$$N_6$$N_7$$N_8$$N_9$
$N_1$210120000
$N_2$131022000
$N_3$011001000
$N_4$100320120
$N_5$220262022
$N_6$021023001
$N_7$000100110
$N_8$000220132
$N_9$000021022

让我逐个验证几个之前可能出错的条目：

$N_2$ 与 $N_6$：

  • $N_2$ 在 {Tri1, Tri3, Tri4}
  • $N_6$ 在 {Tri3, Tri4, Tri7}
  • 交集 = {Tri3, Tri4} → 2个共享单元

之前我写成了0，这是错误的！$M_{2,6} = 2\mu \neq 0$。

$N_4$ 与 $N_8$：

  • $N_4$ 在 {Tri2, Tri5, Tri6}
  • $N_8$ 在 {Tri5, Tri6, Tri8}
  • 交集 = {Tri5, Tri6} → 2个

$M_{4,8} = 2\mu$。之前在刚度矩阵中 $K_{4,8}=0$，但质量矩阵中 $M_{4,8}\neq 0$！

$N_8$ 与 $N_9$：

  • $N_8$ 在 {Tri5, Tri6, Tri8}
  • $N_9$ 在 {Tri7, Tri8}
  • 交集 = {Tri8} → 1个

但是等一下——让我确认 $N_8$ 是否在 Tri7 中。Tri7 的节点是 $(N_5, N_6, N_9)$，不包含 $N_8$。✓

再确认 $N_9$ 是否在 Tri5 或 Tri6 中。Tri5=$(N_4,N_5,N_8)$，Tri6=$(N_4,N_8,N_7)$。都不包含 $N_9$。✓

所以 $M_{8,9}$：交集 = {Tri8} → 1个共享单元 → $M_{8,9} = 1\mu$。

等等，我上面的表写的是2。让我再仔细看。

重新检查 $N_8$ 与 $N_9$：

  • $N_8$ 在 {Tri5, Tri6, Tri8}
  • $N_9$ 在 {Tri7, Tri8}
  • 交集 = {Tri8}
  • 共享1个单元

我表格里写了2，这是错误的。修正为 1。

$N_5$ 与 $N_9$：

  • $N_5$ 在 {Tri1, Tri2, Tri4, Tri5, Tri7, Tri8}
  • $N_9$ 在 {Tri7, Tri8}
  • 交集 = {Tri7, Tri8} → 2个 ✓

$N_4$ 与 $N_7$：

  • $N_4$ 在 {Tri2, Tri5, Tri6}
  • $N_7$ 在 {Tri6}
  • 交集 = {Tri6} → 1个 ✓

$N_7$ 与 $N_8$：

  • $N_7$ 在 {Tri6}
  • $N_8$ 在 {Tri5, Tri6, Tri8}
  • 交集 = {Tri6} → 1个 ✓

$N_6$ 与 $N_9$：

  • $N_6$ 在 {Tri3, Tri4, Tri7}
  • $N_9$ 在 {Tri7, Tri8}
  • 交集 = {Tri7} → 1个 ✓

5.9 修正后的共享单元数表

$(J,I)$$N_1$$N_2$$N_3$$N_4$$N_5$$N_6$$N_7$$N_8$$N_9$
$N_1$210120000
$N_2$131022000
$N_3$011001000
$N_4$100320120
$N_5$220262022
$N_6$021023001
$N_7$000100110
$N_8$000220131
$N_9$000021012

修正项用粗体标注。

5.10 质量矩阵的元素值

质量矩阵元素的规则：

$$ M_{JI} = \text{（共享单元数）} \times \begin{cases} 2\mu & J = I \text{（对角）}\\ \mu & J \neq I \text{（非对角）} \end{cases} $$

其中 $\mu = \frac{\beta A^e}{12} = \frac{(3+j)}{96}$

但这个规则有一个微妙之处需要验证！

对于质量矩阵，在单元 $e$ 中，$M^e_{ji}$ 的值取决于 $j$ 和 $i$ 是否相等（局部编号）。两个不同的全局节点在同一单元中一定有不同的局部编号，所以非对角元素确实用 $\mu$（而非 $2\mu$）。

但对角元素：同一个全局节点在一个单元中只出现一次，对应一个局部编号 $k$，贡献 $M^e_{kk} = 2\mu$。

所以上面的规则是正确的。

5.11 完整修正后的全局质量矩阵

$$ \boxed{ [M]_{9\times 9} = \mu\begin{pmatrix} 4 & 1 & 0 & 1 & 2 & 0 & 0 & 0 & 0 \\ 1 & 6 & 1 & 0 & 2 & 2 & 0 & 0 & 0 \\ 0 & 1 & 2 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 6 & 2 & 0 & 1 & 2 & 0 \\ 2 & 2 & 0 & 2 & 12 & 2 & 0 & 2 & 2 \\ 0 & 2 & 1 & 0 & 2 & 6 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & 2 & 1 & 0 \\ 0 & 0 & 0 & 2 & 2 & 0 & 1 & 6 & 1 \\ 0 & 0 & 0 & 0 & 2 & 1 & 0 & 1 & 4 \end{pmatrix} } $$

验证：

  1. 对称性：$M_{JI} = M_{IJ}$ ✓

  2. 每行和：$\sum_I M_{JI}$ 应与节点的"面积份额"成正比

    • $N_5$（内部）：行和 $= (2+2+0+2+12+2+0+2+2)\mu = 24\mu$
    • $N_1$（角点）：行和 $= (4+1+0+1+2+0+0+0+0)\mu = 8\mu$
    • $N_2$（边中）：行和 $= (1+6+1+0+2+2+0+0+0)\mu = 12\mu$
    • 比例 $N_5:N_2:N_1 = 24:12:8 = 3:1.5:1$

    对于均匀网格，$N_5$（内部）周围面积最大，$N_1$（角点）最小，定性合理 ✓

  3. 总和：$\sum_{J,I}M_{JI} = \beta\sum_e A^e\int_{\Omega^e}\left(\sum_j\lambda_j\right)\left(\sum_i\lambda_i\right)d\Omega = \beta\sum_e A^e\cdot A^e = \beta\cdot A^e\cdot M\cdot A^e$

    不对——应该用 $\sum_{J,I}M_{JI} = \beta\int_\Omega\left(\sum_J\mathcal{N}_J\right)\left(\sum_I\mathcal{N}_I\right)d\Omega = \beta\int_\Omega 1\cdot 1\,d\Omega = \beta\cdot|\Omega| = \beta\cdot 1 = 3+j$

    计算：总和 $= \mu\times(8+12+4+12+24+12+4+12+8) = 96\mu = 96\times\frac{3+j}{96} = 3+j$ ✓✓✓

5.12 全局 Robin 矩阵

$$ [R^b] = \frac{\gamma L^b}{6}\begin{pmatrix}2&1\\1&2\end{pmatrix} = \frac{j\times 0.5}{6}\begin{pmatrix}2&1\\1&2\end{pmatrix} = \nu\begin{pmatrix}2&1\\1&2\end{pmatrix} $$

其中 $\nu = \frac{j}{12}$

逐条边界段散布：

段全局节点贡献位置
$\Gamma^1$: $(N_1,N_2)$$R_{1,1}+=2\nu,\;R_{1,2}+=\nu,\;R_{2,1}+=\nu,\;R_{2,2}+=2\nu$
$\Gamma^2$: $(N_2,N_3)$$R_{2,2}+=2\nu,\;R_{2,3}+=\nu,\;R_{3,2}+=\nu,\;R_{3,3}+=2\nu$
$\Gamma^3$: $(N_3,N_6)$$R_{3,3}+=2\nu,\;R_{3,6}+=\nu,\;R_{6,3}+=\nu,\;R_{6,6}+=2\nu$
$\Gamma^4$: $(N_6,N_9)$$R_{6,6}+=2\nu,\;R_{6,9}+=\nu,\;R_{9,6}+=\nu,\;R_{9,9}+=2\nu$
$\Gamma^5$: $(N_9,N_8)$$R_{9,9}+=2\nu,\;R_{9,8}+=\nu,\;R_{8,9}+=\nu,\;R_{8,8}+=2\nu$
$\Gamma^6$: $(N_8,N_7)$$R_{8,8}+=2\nu,\;R_{8,7}+=\nu,\;R_{7,8}+=\nu,\;R_{7,7}+=2\nu$
$\Gamma^7$: $(N_7,N_4)$$R_{7,7}+=2\nu,\;R_{7,4}+=\nu,\;R_{4,7}+=\nu,\;R_{4,4}+=2\nu$
$\Gamma^8$: $(N_4,N_1)$$R_{4,4}+=2\nu,\;R_{4,1}+=\nu,\;R_{1,4}+=\nu,\;R_{1,1}+=2\nu$

汇总对角元素：

节点所在边界段数$R_{II}$
$N_1$（角点）2（$\Gamma^1, \Gamma^8$）$2\times 2\nu = 4\nu$
$N_2$（边中）2（$\Gamma^1, \Gamma^2$）$2\times 2\nu = 4\nu$
$N_3$（角点）2（$\Gamma^2, \Gamma^3$）$4\nu$
$N_4$（边中）2（$\Gamma^7, \Gamma^8$）$4\nu$
$N_5$（内部）0$0$
$N_6$（边中）2（$\Gamma^3, \Gamma^4$）$4\nu$
$N_7$（角点）2（$\Gamma^6, \Gamma^7$）$4\nu$
$N_8$（边中）2（$\Gamma^5, \Gamma^6$）$4\nu$
$N_9$（角点）2（$\Gamma^5, \Gamma^4$）$4\nu$

所有边界节点的 $R_{II}$ 都等于 $4\nu$，因为每个边界节点恰好在2条边界段上。

汇总非对角元素（只列非零）：

$$ R_{1,2} = \nu, \quad R_{1,4} = \nu $$$$ R_{2,3} = \nu $$$$ R_{3,6} = \nu $$$$ R_{4,7} = \nu $$$$ R_{6,9} = \nu $$$$ R_{7,8} = \nu $$$$ R_{8,9} = \nu $$

加上对称位置。

完整的全局 Robin 矩阵：

$$ \boxed{ [R]_{9\times 9} = \nu\begin{pmatrix} 4 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 1 & 4 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 4 & 0 & 0 & 1 & 0 & 0 & 0 \\ 1 & 0 & 0 & 4 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 4 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & 4 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 4 & 1 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 1 & 4 \end{pmatrix} } $$

验证：

  1. 第5行和第5列全部为零（$N_5$ 是内部节点）✓

  2. 对称 ✓

  3. 总和：$\sum_{J,I}R_{JI} = \gamma\oint_{\Gamma_2}\left(\sum_J\mathcal{N}_J\right)\left(\sum_I\mathcal{N}_I\right)d\Gamma = \gamma\oint_{\Gamma_2}1\,d\Gamma = \gamma\cdot L_{total} = j\times 4 = 4j$

    计算：总和 $= \nu\times(4+1+1+1+4+1+1+4+1+1+1+4+0+1+4+1+1+1+4+1+1+4+1+1+1+4)$

    简便方法：每行和分别为 $6,6,6,6,0,6,6,6,6$，总和 $= 48\nu = 48\times\frac{j}{12} = 4j$ ✓✓✓

六、全局右端项

6.1 体积源项

$$ f^e_j = \frac{f\cdot A^e}{3} = \frac{1\times 0.125}{3} = \frac{1}{24} $$

每个节点收到的体积源 $= \frac{1}{24}\times$（所在单元数）

节点所在单元数$b_J^{(f)}$
$N_1$2$2/24 = 1/12$
$N_2$3$3/24 = 1/8$
$N_3$1$1/24$
$N_4$3$3/24 = 1/8$
$N_5$6$6/24 = 1/4$
$N_6$3$3/24 = 1/8$
$N_7$1$1/24$
$N_8$3$3/24 = 1/8$
$N_9$2$2/24 = 1/12$

验证：$\sum_J b_J^{(f)} = \frac{2+3+1+3+6+3+1+3+2}{24} = \frac{24}{24} = 1 = f\cdot|\Omega|$ ✓

6.2 边界源项

$q = 0$，所以 $b_J^{(q)} = 0$ 对所有 $J$。

6.3 完整右端项

$$ \boxed{ \{b\} = \frac{1}{24}\begin{pmatrix}2\\3\\1\\3\\6\\3\\1\\3\\2\end{pmatrix} } $$

七、最终全局方程组

$$ \boxed{ [A]\{\phi\} = \{b\} } $$$$ [A] = [K] + [M] + [R] $$

7.1 数值代入

$$ \mu = \frac{3+j}{96}, \quad \nu = \frac{j}{12} $$$$ A_{JI} = K_{JI} + \mu\cdot\tilde{M}_{JI} + \nu\cdot\tilde{R}_{JI} $$

其中 $\tilde{M}_{JI}$ 和 $\tilde{R}_{JI}$ 是前面矩阵中的整数系数。

以几个典型元素为例：

$A_{5,5}$（内部节点对角元）：

$$ A_{5,5} = 6.0 + 12\mu + 0 = 6.0 + \frac{12(3+j)}{96} + 0 = 6.0 + \frac{3+j}{8} $$$$ = 6.0 + 0.375 + 0.125j = 6.375 + 0.125j $$

$A_{1,1}$（角点对角元）：

$$ A_{1,1} = 1.5 + 4\mu + 4\nu = 1.5 + \frac{4(3+j)}{96} + \frac{4j}{12} $$$$ = 1.5 + \frac{3+j}{24} + \frac{j}{3} = 1.5 + 0.125 + \frac{j}{24} + \frac{j}{3} $$$$ = 1.625 + \frac{j+8j}{24} = 1.625 + \frac{9j}{24} = 1.625 + 0.375j $$

$A_{5,2}$（内部-边界耦合，非对角）：

$$ A_{5,2} = -2.0 + 2\mu + 0 = -2.0 + \frac{2(3+j)}{96} = -2.0 + \frac{3+j}{48} $$$$ = -2.0 + 0.0625 + 0.02083j = -1.9375 + 0.02083j $$

$A_{2,6}$（两个边界节点之间）：

$$ A_{2,6} = K_{2,6} + M_{2,6}\cdot\mu + R_{2,6}\cdot\nu $$

$K_{2,6} = 0$（从刚度矩阵可以看到），$M_{2,6} = 2\mu$，$R_{2,6} = 0$（$N_2$ 和 $N_6$ 不共享边界段）

$$ A_{2,6} = 0 + 2\mu + 0 = \frac{2(3+j)}{96} = \frac{3+j}{48} $$

虽然刚度矩阵中 $K_{2,6} = 0$，但质量矩阵使得 $A_{2,6} \neq 0$！

7.2 方程组的结构特征

$$ \text{方程组大小：} 9\times 9 $$$$ \text{矩阵性质：复数、对称（因为 } \bar{\bar{\alpha}}\text{ 对角且 } \beta,\gamma\text{ 是标量）、稀疏} $$

非零元素的模式：

      N1   N2   N3   N4   N5   N6   N7   N8   N9
N1  [  •    •    ·    •    •    ·    ·    ·    ·  ]
N2  [  •    •    •    ·    •    •    ·    ·    ·  ]
N3  [  ·    •    •    ·    ·    •    ·    ·    ·  ]
N4  [  •    ·    ·    •    •    ·    •    •    ·  ]
N5  [  •    •    ·    •    •    •    ·    •    •  ]
N6  [  ·    •    •    ·    •    •    ·    ·    •  ]
N7  [  ·    ·    ·    •    ·    ·    •    •    ·  ]
N8  [  ·    ·    ·    •    •    ·    •    •    •  ]
N9  [  ·    ·    ·    ·    •    •    ·    •    •  ]

注意：$A_{2,6}\neq 0$（修正后），$A_{4,8}\neq 0$，但 $K_{2,6}=K_{4,8}=0$。这些位置在刚度矩阵中为零，但在质量矩阵中非零。

这意味着全局矩阵 $[A]$ 的稀疏模式是三个子矩阵稀疏模式的并集，比任何单个子矩阵更"稠密"。

八、内部节点 vs 边界节点——方程的物理含义

8.1 内部节点 $N_5$ 的方程（第5行）

$$ A_{5,1}\phi_1 + A_{5,2}\phi_2 + A_{5,4}\phi_4 + A_{5,5}\phi_5 + A_{5,6}\phi_6 + A_{5,8}\phi_8 + A_{5,9}\phi_9 = b_5 $$

特征：

  • 没有 Robin 边界贡献（$R_{5,\cdot}$ 全为零）
  • 只有体积刚度 + 体积质量
  • 与所有6个相邻节点耦合
  • 对角元素最大（$6.375+0.125j$）

8.2 角点 $N_1$ 的方程（第1行）

$$ A_{1,1}\phi_1 + A_{1,2}\phi_2 + A_{1,4}\phi_4 + A_{1,5}\phi_5 = b_1 $$

特征：

  • 有 Robin 边界贡献（$R_{1,1}, R_{1,2}, R_{1,4}$ 非零）
  • 只与3个节点耦合（$N_2, N_4, N_5$）
  • Robin 条件已经"自然地"嵌入了方程，不需要额外处理

8.3 边中点 $N_2$ 的方程（第2行）

$$ A_{2,1}\phi_1 + A_{2,2}\phi_2 + A_{2,3}\phi_3 + A_{2,5}\phi_5 + A_{2,6}\phi_6 = b_2 $$

特征：

  • Robin 贡献在 $A_{2,1}, A_{2,2}, A_{2,3}$ 中（来自 $\Gamma^1, \Gamma^2$）
  • 与 $N_6$ 的耦合仅来自质量矩阵，不来自刚度和Robin

8.4 与三维的完整对比

特征2D节点元（本文）3D棱边元
DOF位置节点上棱边上
DOF物理含义$\phi_I = \phi(x_I,y_I)$，节点处的场值$E_I = \int_{E_I}\mathbf{E}\cdot d\mathbf{l}$，棱边切向线积分
每个单元DOF数3（三角形）6（四面体）
基函数类型标量 $N_k^{(e)} = \lambda_k$矢量 $\mathbf{N}_k^{(e)} = \lambda_a\nabla\lambda_b - \lambda_b\nabla\lambda_a$
基函数连续性场本身连续（$C^0$）切向连续，法向可以不连续
方向符号不需要必须有 $s_k^{(e)} = \pm 1$
全局-局部关系$\mathcal{N}_I\big\|_{\Omega^e} = N_{i}^{(e)}$$\boldsymbol{\mathcal{N}}_I\big\|_{\Omega^e} = s_i^{(e)}\mathbf{N}_i^{(e)}$
组装公式$A_{JI} \mathrel{+}= K^e_{ji}$$A_{JI} \mathrel{+}= s_j^{(e)}s_i^{(e)}K^e_{ji}$
关联矩阵$C^{(e)}_{I,k} = \delta_{I,\text{map}(e,k)}$$C^{(e)}_{I,k} = s_k^{(e)}\delta_{I,\text{map}(e,k)}$
全局矩阵公式$[A]=\sum_e[L^{(e)}]^T[K^e][L^{(e)}]$$[A]=\sum_e[P^{(e)}]^T[K^e][P^{(e)}]$
$L^{(e)}$ 是纯布尔矩阵$P^{(e)} = \Sigma^{(e)}L^{(e)}$ 含符号

刚度矩阵的类比：

2D节点元3D棱边元
公式$K^e_{ji} = \frac{1}{4A^e}\mathbf{g}_j^T\bar{\bar{D}}\,\mathbf{g}_i$$S^e_{ji} = 4(\nabla\lambda_{a_j}\times\nabla\lambda_{b_j})\cdot(\nabla\lambda_{a_i}\times\nabla\lambda_{b_i})V^e$
$\mathbf{g}_k = (b_k, c_k)^T$$\mathbf{c}_k = \nabla\lambda_{a_k}\times\nabla\lambda_{b_k}$
本质梯度的加权内积 × 面积旋度的内积 × 体积
性质常数被积函数 → 精确积分常数被积函数 → 精确积分
各向异性通过 $\bar{\bar{D}}=\text{diag}(\alpha_x,\alpha_y)$通过 $D_{pq}=(\nabla\lambda_p)^T\bar{\bar{\varepsilon}}_r(\nabla\lambda_q)$

质量矩阵的类比：

2D节点元3D棱边元
公式$M^e_{ji} = \beta\int_{\Omega^e}\lambda_j\lambda_i\,d\Omega$$T^e_{ji} = \int_{\Omega^e}\mathbf{N}_j\cdot\bar{\bar{\varepsilon}}_r\mathbf{N}_i\,dV$
积分公式$\int\lambda_j\lambda_i = \frac{(1+\delta_{ji})A^e}{12}$四项展开，用 $I_{pq}=\frac{(1+\delta_{pq})V^e}{20}$
对角/非对角比$2:1$依赖于具体节点对
各向异性标量 $\beta$ 直接提出积分张量 $\bar{\bar{\varepsilon}}_r$ 产生复杂的四项展开

边界条件的类比：

2D节点元（Robin）3D棱边元（ABC）
条件$(\bar{\bar{\alpha}}\nabla\phi)\cdot\hat{n}+\gamma\phi=q$$\hat{n}\times(\nabla\times\mathbf{E})=jk_0\hat{n}\times(\hat{n}\times\mathbf{E})$
类型自然边界条件（代入弱形式）自然边界条件（代入弱形式）
产生的矩阵$R^b_{ji}=\frac{\gamma L^b}{6}\begin{pmatrix}2&1\\1&2\end{pmatrix}$$\tilde{B}^f_{\beta\alpha}$（面上3×3矩阵）
积分域1D线段（2个节点）2D三角形面（3条棱边）
积分公式$\frac{(1+\delta)L^b}{6}$$\frac{(1+\delta)A^f}{12}$
仅影响边界节点（内部节点无贡献）边界棱边（内部棱边无贡献）

维度递减规律：

            体积积分          面/边界积分
            ─────────         ──────────
3D四面体:   ∫λ_iλ_j dV        ∫λ_iλ_j dS
            (1+δ)V^e/20      (1+δ)A^f/12

2D三角形:   ∫λ_iλ_j dΩ       ∫N_iN_j dΓ
            (1+δ)A^e/12      (1+δ)L^b/6

1D线段:     ∫N_iN_j dx        N_iN_j|端点
            (1+δ)L^e/6       δ_{ij}

规律:  d维:  (1+δ)·|Ω|/((d+1)(d+2)/2·...)
分母序列:    20, 12, 6 → 即 (d+2)!/2 的某种形式

支撑域筛选的类比：

2D节点元3D棱边元
支撑域包含节点 $I$ 的所有三角形包含棱边 $I$ 的所有四面体
典型共享单元数5-7个（内部节点）5-7个（内部棱边）
$A_{JI}\neq 0$ 条件$I$ 和 $J$ 共享至少一个三角形$I$ 和 $J$ 共享至少一个四面体
每行非零元素~7个~15-20个
零的原因$\mathcal{N}_I\big\|_{\Omega^e}=0$ 若 $I\notin e$$\boldsymbol{\mathcal{N}}_I\big\|_{\Omega^e}=\mathbf{0}$ 若 $I\notin e$
边界额外筛选Robin仅涉及边界线段上的节点ABC仅涉及面上的棱边
右端额外筛选源函数为零的区域不贡献$(\bar{\bar{\varepsilon}}_r-\bar{\bar{I}})=0$ 的区域不贡献

组装过程的完全类比：

2D节点元：                          3D棱边元：

(1) 定义全局节点                    (1) 定义全局棱边
    I = 1,...,N                         I = 1,...,N_E
    
(2) 建立映射表                      (2) 建立映射表 + 符号表
    map(e,k) → I                        map(e,k) → I
    无符号                               sign(e,k) → ±1

(3) 全局展开                        (3) 全局展开
    φ_h = Σ φ_I N_I                     E_h = Σ E_I N_I

(4) 代入弱形式                      (4) 代入弱形式
    Σ A_{JI} φ_I = b_J                  Σ A_{JI} E_I = b_J

(5) 全域积分→单元积分               (5) 全域积分→单元积分
    A_{JI} = Σ_e (...)                  A_{JI} = Σ_e (...)

(6) 支撑域筛选                      (6) 支撑域筛选
    N_I|_{Ω^e} = 0 若 I∉e              N_I|_{Ω^e} = 0 若 I∉e

(7) 全局→局部                       (7) 全局→局部
    N_I|_{Ω^e} = N_i^(e)               N_I|_{Ω^e} = s_i^(e) N_i^(e)
    无符号!                              有符号!

(8) 局部矩阵计算                    (8) 局部矩阵计算
    K^e_{ji}, M^e_{ji}                  S^e_{ji}, T^e_{ji}
    纯局部，不涉及全局                   纯局部，不涉及全局

(9) 局部→全局（组装）               (9) 局部→全局（组装）
    A_{JI} += K^e_{ji}                  A_{JI} += s_j s_i K^e_{ji}
    无符号!                              有符号!

(10) 边界贡献                       (10) 边界贡献
     Robin: 线段上2×2               ABC: 三角形上3×3
     A_{JI} += R^b_{ji}             A_{JI} += s_α s_β B^f_{βα}
     无符号!                          有符号!

(11) 求解 [A]{φ} = {b}             (11) 求解 [A]{E} = {b}

错误风险的对比：

错误类型2D节点元3D棱边元
忘记方向符号不可能（没有符号）致命错误，解完全错误
映射表错误会出错会出错
边界条件遗漏Robin项丢失ABC项丢失
积分公式系数错会出错会出错
伪解（spurious modes）不会产生节点元会产生，棱边元不会

总结：2D节点元是3D棱边元的"简化版"。两者在全局→局部→全局的逻辑框架上完全一致，区别仅在于：

  1. 棱边元多了方向符号 $s_k^{(e)}$
  2. 棱边元的基函数是矢量而非标量
  3. 棱边元的质量矩阵因各向异性张量而更复杂
  4. 棱边元保证了矢量场的正确物理行为（无伪解）

九、完整的代码实现

import numpy as np

# ═══════════════════════════════════════
# 参数
# ═══════════════════════════════════════
alpha_x = 1.0
alpha_y = 2.0
beta = 3.0 + 1j
gamma = 1j
f_source = 1.0
q_source = 0.0

# ═══════════════════════════════════════
# 网格数据
# ═══════════════════════════════════════
node_coords = np.array([
    [0.0, 0.0],   # N1
    [0.5, 0.0],   # N2
    [1.0, 0.0],   # N3
    [0.0, 0.5],   # N4
    [0.5, 0.5],   # N5
    [1.0, 0.5],   # N6
    [0.0, 1.0],   # N7
    [0.5, 1.0],   # N8
    [1.0, 1.0],   # N9
])

# 单元连接表（0-based）
elem_nodes = np.array([
    [0, 1, 4],  # Tri1: N1,N2,N5
    [0, 4, 3],  # Tri2: N1,N5,N4
    [1, 2, 5],  # Tri3: N2,N3,N6
    [1, 5, 4],  # Tri4: N2,N6,N5
    [3, 4, 7],  # Tri5: N4,N5,N8
    [3, 7, 6],  # Tri6: N4,N8,N7
    [4, 5, 8],  # Tri7: N5,N6,N9
    [4, 8, 7],  # Tri8: N5,N9,N8
])

# 边界线段（0-based）
boundary_segs = np.array([
    [0, 1],  # Γ1: N1-N2
    [1, 2],  # Γ2: N2-N3
    [2, 5],  # Γ3: N3-N6
    [5, 8],  # Γ4: N6-N9
    [8, 7],  # Γ5: N9-N8
    [7, 6],  # Γ6: N8-N7
    [6, 3],  # Γ7: N7-N4
    [3, 0],  # Γ8: N4-N1
])

N_nodes = len(node_coords)
N_elem = len(elem_nodes)
N_bseg = len(boundary_segs)

# ═══════════════════════════════════════
# 全局矩阵初始化
# ═══════════════════════════════════════
A = np.zeros((N_nodes, N_nodes), dtype=complex)
b = np.zeros(N_nodes, dtype=complex)

# ═══════════════════════════════════════
# 体积单元循环
# ═══════════════════════════════════════
for e in range(N_elem):
    # 获取3个节点坐标
    n = elem_nodes[e]
    x1, y1 = node_coords[n[0]]
    x2, y2 = node_coords[n[1]]
    x3, y3 = node_coords[n[2]]
    
    # 面积
    two_A = (x1 - x3)*(y2 - y3) - (x2 - x3)*(y1 - y3)
    Ae = abs(two_A) / 2.0
    
    # b_k, c_k 系数
    b_coeff = np.array([y2 - y3, y3 - y1, y1 - y2])
    c_coeff = np.array([x3 - x2, x1 - x3, x2 - x1])
    
    # 局部刚度矩阵 3×3
    Ke = np.zeros((3, 3))
    for j in range(3):
        for i in range(3):
            Ke[j, i] = (alpha_x * b_coeff[j] * b_coeff[i] +
                        alpha_y * c_coeff[j] * c_coeff[i]) / (4.0 * Ae)
    
    # 局部质量矩阵 3×3
    Me = np.zeros((3, 3), dtype=complex)
    for j in range(3):
        for i in range(3):
            if j == i:
                Me[j, i] = beta * Ae / 6.0
            else:
                Me[j, i] = beta * Ae / 12.0
    
    # 局部右端项
    fe = f_source * Ae / 3.0 * np.ones(3)
    
    # 组装到全局（无符号！）
    for j in range(3):
        J = n[j]  # 全局行号
        b[J] += fe[j]
        for i in range(3):
            I = n[i]  # 全局列号
            A[J, I] += Ke[j, i] + Me[j, i]

# ═══════════════════════════════════════
# 边界线段循环
# ═══════════════════════════════════════
for seg in range(N_bseg):
    n = boundary_segs[seg]
    x1, y1 = node_coords[n[0]]
    x2, y2 = node_coords[n[1]]
    
    Lb = np.sqrt((x2-x1)**2 + (y2-y1)**2)
    
    # 局部 Robin 矩阵 2×2
    Rb = gamma * Lb / 6.0 * np.array([[2, 1], [1, 2]])
    
    # 局部边界源 2×1
    qb = q_source * Lb / 2.0 * np.ones(2)
    
    # 组装到全局（无符号！）
    for j in range(2):
        J = n[j]
        b[J] += qb[j]
        for i in range(2):
            I = n[i]
            A[J, I] += Rb[j, i]

# ═══════════════════════════════════════
# 求解
# ═══════════════════════════════════════
phi = np.linalg.solve(A, b)

print("解 φ：")
for i in range(N_nodes):
    print(f"  N{i+1} ({node_coords[i][0]:.1f}, {node_coords[i][1]:.1f}): "
          f"φ = {phi[i]:.6f}")

# ═══════════════════════════════════════
# 验证
# ═══════════════════════════════════════
residual = np.linalg.norm(A @ phi - b)
print(f"\n残差: {residual:.2e}")
print(f"矩阵对称性误差: {np.linalg.norm(A - A.T):.2e}")

# 验证刚度矩阵行和为零
K_only = np.zeros((N_nodes, N_nodes))
for e in range(N_elem):
    n = elem_nodes[e]
    x1, y1 = node_coords[n[0]]
    x2, y2 = node_coords[n[1]]
    x3, y3 = node_coords[n[2]]
    two_A = (x1-x3)*(y2-y3) - (x2-x3)*(y1-y3)
    Ae = abs(two_A)/2.0
    bc = np.array([y2-y3, y3-y1, y1-y2])
    cc = np.array([x3-x2, x1-x3, x2-x1])
    for j in range(3):
        for i in range(3):
            K_only[n[j], n[i]] += (alpha_x*bc[j]*bc[i] + 
                                    alpha_y*cc[j]*cc[i])/(4*Ae)

print(f"刚度矩阵最大行和: {np.max(np.abs(np.sum(K_only, axis=1))):.2e}")