这篇笔记在做什么
这一篇对应原始会话里最完整的一条三维电磁有限元推导链。我保留了四面体一阶 Nedelec 棱边元、$e^{-j\omega t}$ 约定、ABC 面矩阵、全局与局部映射、质量矩阵和代码框架的细致展开。
如果你关心的是“从 Maxwell 方程一路走到可实现的单元矩阵和组装算法”，这一篇就是主干正文。
三维各向异性介质散射问题有限元方法
——四面体棱边元 · ABC边界 · $e^{-j\omega t}$约定 · 完整推导
一、物理问题与几何设置
各向异性无磁电介质立方体，被平面波照射：
E^inc (平面波)
↓ ↓ ↓
╔══════════════════════════╗
║ Ω₀ (自由空间, ε_r=I) ║ ← S_ABC (外边界, 球/长方体)
║ ║
║ ┌──────────┐ ║
║ │ Ω_d │ ║
║ │ 各向异性 │ ║
║ │ 介质体 │ ║
║ │ ε̄̄_r(张量)│ ║
║ └──────────┘ ║
║ S_d (介质面) ║
║ ║
╚══════════════════════════╝
μ_r = 1 (全空间无磁)
材料参数——各向异性介电张量：
$$ \bar{\bar{\varepsilon}}_r = \begin{pmatrix} \varepsilon_{xx} & \varepsilon_{xy} & \varepsilon_{xz} \\ \varepsilon_{yx} & \varepsilon_{yy} & \varepsilon_{yz} \\ \varepsilon_{zx} & \varepsilon_{zy} & \varepsilon_{zz} \end{pmatrix} $$区域划分：
$$ \bar{\bar{\varepsilon}}_r(\mathbf{r}) = \begin{cases} \bar{\bar{\varepsilon}}_r^{(d)} & \mathbf{r} \in \Omega_d \text{（介质内部）} \\ \bar{\bar{I}} & \mathbf{r} \in \Omega_0 \text{（自由空间）} \end{cases} $$二、Maxwell方程与波动方程（$e^{-j\omega t}$）
2.1 时谐Maxwell方程
采用 $e^{-j\omega t}$ 因子约定：
$$ \nabla \times \mathbf{E} = j\omega\mu_0 \mathbf{H} $$$$ \nabla \times \mathbf{H} = -j\omega \bar{\bar{\varepsilon}}_0 \bar{\bar{\varepsilon}}_r \mathbf{E} $$注意符号：$e^{-j\omega t}$ 与 $e^{+j\omega t}$ 的区别仅在于 $j \leftrightarrow -j$。
2.2 矢量波动方程
由 $\mathbf{H} = \frac{1}{j\omega\mu_0}\nabla\times\mathbf{E}$，代入第二式（$\mu_r=1$）：
$$ \nabla \times \left(\nabla \times \mathbf{E}\right) - k_0^2 \bar{\bar{\varepsilon}}_r \mathbf{E} = \mathbf{0} $$其中 $k_0 = \omega\sqrt{\mu_0\varepsilon_0}$。
2.3 散射场公式
设 $\mathbf{E}^{tot} = \mathbf{E}^{inc} + \mathbf{E}^{sca}$。
入射场满足自由空间方程：
$$ \nabla\times(\nabla\times\mathbf{E}^{inc}) - k_0^2\mathbf{E}^{inc} = \mathbf{0} $$总场满足：
$$ \nabla\times(\nabla\times\mathbf{E}^{tot}) - k_0^2 \bar{\bar{\varepsilon}}_r \mathbf{E}^{tot} = \mathbf{0} $$代入 $\mathbf{E}^{tot} = \mathbf{E}^{inc} + \mathbf{E}^{sca}$：
$$ \nabla\times(\nabla\times\mathbf{E}^{sca}) - k_0^2 \bar{\bar{\varepsilon}}_r \mathbf{E}^{sca} = k_0^2(\bar{\bar{\varepsilon}}_r - \bar{\bar{I}})\mathbf{E}^{inc} $$关键：右端项只在介质内部 $\Omega_d$ 非零（因为自由空间中 $\bar{\bar{\varepsilon}}_r = \bar{\bar{I}}$）。
定义：
$$ \boxed{ \nabla\times(\nabla\times\mathbf{E}^{sca}) - k_0^2 \bar{\bar{\varepsilon}}_r \mathbf{E}^{sca} = \mathbf{f} } $$$$ \mathbf{f}(\mathbf{r}) = \begin{cases} k_0^2(\bar{\bar{\varepsilon}}_r - \bar{\bar{I}})\mathbf{E}^{inc} & \mathbf{r} \in \Omega_d \\ \mathbf{0} & \mathbf{r} \in \Omega_0 \end{cases} $$三、一阶ABC边界条件（$e^{-j\omega t}$）
3.1 推导
远场散射波近似为局部平面波，传播方向 $\hat{r} \approx \hat{n}$（外法线）：
$$ \mathbf{E}^{sca} \sim \frac{e^{jk_0 r}}{r}\mathbf{F}(\theta,\phi) $$注意：$e^{-j\omega t}$ 下，外行波空间因子为 $e^{+jk_0 r}$。
远场关系：
$$ \nabla\times\mathbf{E}^{sca} \approx jk_0 \hat{n}\times\mathbf{E}^{sca} \quad \text{（在 } S_{ABC} \text{ 上）} $$3.2 Silver-Müller ABC
$$ \boxed{ \hat{n}\times(\nabla\times\mathbf{E}^{sca}) = jk_0\,\hat{n}\times(\hat{n}\times\mathbf{E}^{sca}) \quad \text{在 } S_{ABC} \text{ 上} } $$对比：若用 $e^{+j\omega t}$，则符号取负：$\hat{n}\times(\nabla\times\mathbf{E}^{sca}) = -jk_0\,\hat{n}\times(\hat{n}\times\mathbf{E}^{sca})$。
四、弱形式（变分形式）完整推导
4.1 加权残量
用矢量测试函数 $\mathbf{T}$ 对波动方程做内积：
$$ \int_\Omega \mathbf{T}\cdot\left[\nabla\times(\nabla\times\mathbf{E}^{sca}) - k_0^2\bar{\bar{\varepsilon}}_r\mathbf{E}^{sca}\right]dV = \int_\Omega \mathbf{T}\cdot\mathbf{f}\,dV $$4.2 分部积分（矢量Green恒等式）
利用恒等式：
$$ \nabla\cdot(\mathbf{T}\times\nabla\times\mathbf{E}) = (\nabla\times\mathbf{T})\cdot(\nabla\times\mathbf{E}) - \mathbf{T}\cdot(\nabla\times\nabla\times\mathbf{E}) $$得：
$$ \mathbf{T}\cdot(\nabla\times\nabla\times\mathbf{E}) = (\nabla\times\mathbf{T})\cdot(\nabla\times\mathbf{E}) - \nabla\cdot(\mathbf{T}\times\nabla\times\mathbf{E}) $$对体积积分，应用散度定理：
$$ \int_\Omega \mathbf{T}\cdot(\nabla\times\nabla\times\mathbf{E})\,dV = \int_\Omega (\nabla\times\mathbf{T})\cdot(\nabla\times\mathbf{E})\,dV - \oint_S \hat{n}\cdot(\mathbf{T}\times\nabla\times\mathbf{E})\,dS $$利用标量三重积恒等式：
$$ \hat{n}\cdot(\mathbf{T}\times\nabla\times\mathbf{E}) = (\hat{n}\times\mathbf{T})\cdot(\nabla\times\mathbf{E}) $$因此：
$$ \int_\Omega \mathbf{T}\cdot(\nabla\times\nabla\times\mathbf{E})\,dV = \int_\Omega (\nabla\times\mathbf{T})\cdot(\nabla\times\mathbf{E})\,dV - \oint_S (\hat{n}\times\mathbf{T})\cdot(\nabla\times\mathbf{E})\,dS $$4.3 代入边界条件
在 $S_{ABC}$ 上，由ABC：$\nabla\times\mathbf{E}^{sca} = jk_0\hat{n}\times\mathbf{E}^{sca}$
面积分项变为：
$$ -\oint_{S_{ABC}} (\hat{n}\times\mathbf{T})\cdot(jk_0\hat{n}\times\mathbf{E}^{sca})\,dS = -jk_0 \oint_{S_{ABC}} (\hat{n}\times\mathbf{T})\cdot(\hat{n}\times\mathbf{E}^{sca})\,dS $$四、弱形式（续）
4.4 最终弱形式
$$ \boxed{ \begin{aligned} &\int_\Omega (\nabla\times\mathbf{T})\cdot(\nabla\times\mathbf{E})\,dV \\ &\quad - k_0^2 \int_\Omega \mathbf{T}\cdot\bar{\bar{\varepsilon}}_r\mathbf{E}\,dV \\ &\quad - jk_0 \oint_{S_{ABC}} (\hat{n}\times\mathbf{T})\cdot(\hat{n}\times\mathbf{E})\,dS \\ &= \int_{\Omega_d} \mathbf{T}\cdot k_0^2(\bar{\bar{\varepsilon}}_r - \bar{\bar{I}})\mathbf{E}^{inc}\,dV \end{aligned} } $$各项的物理含义：
| 项 | 表达式 | 含义 |
|---|---|---|
| $[S]$ | $\int_\Omega (\nabla\times\mathbf{T})\cdot(\nabla\times\mathbf{E})\,dV$ | 旋度-旋度项（刚度矩阵） |
| $[T]$ | $-k_0^2\int_\Omega \mathbf{T}\cdot\bar{\bar{\varepsilon}}_r\mathbf{E}\,dV$ | 质量矩阵（含各向异性张量） |
| $[B]$ | $-jk_0\oint_{S_{ABC}}(\hat{n}\times\mathbf{T})\cdot(\hat{n}\times\mathbf{E})\,dS$ | ABC贡献的边界矩阵 |
| $\{f\}$ | $k_0^2\int_{\Omega_d}\mathbf{T}\cdot(\bar{\bar{\varepsilon}}_r-\bar{\bar{I}})\mathbf{E}^{inc}\,dV$ | 入射场激励（仅介质内非零） |
五、四面体单元与棱边元
5.1 四面体单元定义
一个四面体有 4个节点、6条棱边、4个面：
节点4 (n4)
/|\
/ | \
/ | \
e6 / |e5 \ e4
/ | \
/ e3 | \
n1 /______|______\ n3
\ | /
\ | /
e1 \ | / e2
\ | /
\ | /
\ | /
\|/
节点2 (n2)
节点编号： $n_1, n_2, n_3, n_4$
棱边定义与编号约定（起点编号 < 终点编号）：
| 局部棱边号 | 起点 → 终点 | 记法 |
|---|---|---|
| $e_1$ | $n_1 \to n_2$ | (1,2) |
| $e_2$ | $n_1 \to n_3$ | (1,3) |
| $e_3$ | $n_1 \to n_4$ | (1,4) |
| $e_4$ | $n_2 \to n_3$ | (2,3) |
| $e_5$ | $n_2 \to n_4$ | (2,4) |
| $e_6$ | $n_3 \to n_4$ | (3,4) |
每个四面体有6条棱边 → 6个自由度
5.2 体积坐标（重心坐标）
四面体内任意一点用4个体积坐标 $(\lambda_1, \lambda_2, \lambda_3, \lambda_4)$ 描述：
$$ \lambda_1 + \lambda_2 + \lambda_3 + \lambda_4 = 1, \quad \lambda_i \geq 0 $$物理坐标与体积坐标的关系：
$$ \mathbf{r} = \lambda_1 \mathbf{r}_1 + \lambda_2 \mathbf{r}_2 + \lambda_3 \mathbf{r}_3 + \lambda_4 \mathbf{r}_4 $$其中 $\mathbf{r}_i = (x_i, y_i, z_i)^T$ 是第 $i$ 个节点坐标。
体积坐标的梯度（非常重要）：
$$ \nabla\lambda_i = \frac{1}{6V^e}\,\hat{n}_i\,A_i $$更明确地写成分量形式。定义 Jacobian 矩阵：
$$ [J] = \begin{pmatrix} x_1 - x_4 & y_1 - y_4 & z_1 - z_4 \\ x_2 - x_4 & y_2 - y_4 & z_2 - z_4 \\ x_3 - x_4 & y_3 - y_4 & z_3 - z_4 \end{pmatrix} $$$$ V^e = \frac{1}{6}|\det[J]| $$体积坐标梯度由 $[J]$ 的逆给出：
$$ \begin{pmatrix} \nabla\lambda_1 \\ \nabla\lambda_2 \\ \nabla\lambda_3 \end{pmatrix} = [J]^{-T} \begin{pmatrix} \hat{e}_1 \\ \hat{e}_2 \\ \hat{e}_3 \end{pmatrix} $$$$ \nabla\lambda_4 = -\nabla\lambda_1 - \nabla\lambda_2 - \nabla\lambda_3 $$具体展开：设 $[J]^{-1} = (c_{ij})$，则
$$ \nabla\lambda_i = \begin{pmatrix} c_{i1} \\ c_{i2} \\ c_{i3} \end{pmatrix} , \quad i=1,2,3 $$5.3 一阶Nédélec棱边基函数
连接节点 $n_i$ 和 $n_j$（$i 6个基函数完整列表： 关键性质验证： 性质1：沿自身棱边的线积分 = 1 在棱边 $e_{ij}$ 上：$\lambda_i = 1-t$，$\lambda_j = t$，其余 $\lambda_k = 0$，$d\mathbf{l} = (\mathbf{r}_j-\mathbf{r}_i)dt$ 利用 $\nabla\lambda_k\cdot(\mathbf{r}_j-\mathbf{r}_i) = \delta_{kj}-\delta_{ki}$： 性质2：沿其他棱边的线积分 = 0 ✓ 性质3：切向分量跨单元连续 ✓ 利用 $\nabla\times(\phi\nabla\psi) = \nabla\phi\times\nabla\psi$（因为 $\nabla\times(\nabla\psi) = 0$）： 这是一个常矢量（在线性四面体中），大大简化积分！ 完整列出： 单元内电场展开： 代入弱形式，得到 $6\times 6$ 的单元矩阵方程： 由于 $\nabla\times\mathbf{N}_i$ 为常矢量： 其中棱边 $i$ 对应节点对 $(a,b)$，棱边 $j$ 对应节点对 $(c,d)$。 展开成分量： 设 $\nabla\lambda_i = (\alpha_i, \beta_i, \gamma_i)^T$，则 完整的 $6\times6$ 矩阵（每个元素都可以显式写出）： 这是各向异性的核心项。展开基函数： 展开四项： 由于 $\nabla\lambda$ 在线性四面体中是常数，可提出积分： 为保证长公式稳定渲染，先记 展开回原始记号就是： 其中需要计算体积坐标乘积的积分： 四面体中体积坐标积分公式： 因此： 写成紧凑形式： 各向异性的具体展开（以 $T^e_{11}$ 为例）： 棱边1对应 $(a,b)=(1,2)$，$(c,d)=(1,2)$： 即 代入 $I_{11}=I_{22}=V^e/10$，$I_{12}=I_{21}=V^e/20$： 再展开回原始记号： 通用的标量 $(\nabla\lambda_p)^T \bar{\bar{\varepsilon}}_r (\nabla\lambda_q)$ 计算： 设 $\nabla\lambda_p = (\alpha_p, \beta_p, \gamma_p)^T$： 这就是各向异性在FEM中发挥作用的地方。若 $\bar{\bar{\varepsilon}}_r = \varepsilon_r\bar{\bar{I}}$（各向同性），则退化为 $\varepsilon_r\,\nabla\lambda_p\cdot\nabla\lambda_q$。 仅在介质单元内非零。 设入射平面波（$e^{-j\omega t}$ 约定，沿 $\hat{z}$ 传播，$\hat{x}$ 极化）： 注意：$e^{-j\omega t}$ 下空间因子为 $e^{+jk_0 z}$。 令 $\mathbf{g} = (\bar{\bar{\varepsilon}}_r - \bar{\bar{I}})\,E_0\hat{x}$（常矢量，对给定单元），令 $\nabla\lambda_b\cdot\mathbf{g} = g_b$（常数），$\nabla\lambda_a\cdot\mathbf{g} = g_a$（常数）： 对于电尺寸较小的单元（$k_0 h \ll 1$），$e^{jk_0 z} \approx$ 常数 $\approx e^{jk_0 z_c}$（$z_c$ 为单元中心），则： 更精确地，用解析积分或高阶高斯积分。 仅涉及ABC面上的三角形面元 $\Delta^f$。 一个三角形面有3个节点和3条棱边。假设面 $f$ 的三个节点为 $n_p, n_q, n_r$（四面体的第4个节点 $n_s$ 不在面上，即 $\lambda_s = 0$）。 面上 $\lambda_s = 0$，有效基函数为： 其余3条棱边（含节点 $n_s$）在面上切向分量为零，不贡献。 面上的切向基函数： $\hat{n}\times\mathbf{N}_{ij}$ 用二维面积坐标表示。 利用恒等式：$(\hat{n}\times\mathbf{A})\cdot(\hat{n}\times\mathbf{B}) = \mathbf{A}\cdot\mathbf{B} - (\mathbf{A}\cdot\hat{n})(\mathbf{B}\cdot\hat{n})$ 对于棱边元基函数，在面上已经是切向的（$\mathbf{N}_i\cdot\hat{n}=0$ 在面上），所以： 因此： 这与质量矩阵类似，但是面积分。用二维积分公式： $3\times 3$ 面矩阵按照与体积质量矩阵相同的展开方法计算，只是用面积分公式替代体积积分公式。 用5个四面体剖分一个单位立方体，来完整展示编号系统。 8个全局节点： 验证拓扑一致性： 规则：每条棱边用两个全局节点表示，且小编号在前。 遍历所有单元，收集所有不重复的棱边： 验证： 立方体 8节点、5四面体：$E = V - F + C + ...$ Euler公式验证。
或直接数：立方体12条边 + 对角线引入的内部边6条 = 18条。✓ 为什么需要符号？ 全局棱边有一个约定方向（小节点→大节点），但局部棱边也有方向（$n_i \to n_j$, $i 规则： 对于局部棱边 $e_k$ 连接局部节点 $(n_a, n_b)$，其全局节点为 $(N_A, N_B)$： 例子：Tet 2的局部棱边 $e_4$ 符号在组装中的作用： 即，组装时： 其中 $s_i, s_j$ 是对应的符号，$I(i), I(j)$ 是全局棱边号。 我们把最关键的数据结构列成完整表格： 全局方程为 $18\times 18$ 系统 $[A]\{E\} = \{b\}$： 物理原因： 自由度 $E_I$ 代表电场沿全局棱边方向的线积分。如果局部方向与全局方向相反，则： 单元方程： 代入全局： 两边乘以 $s_i$： 所以全局贡献是 $s_i s_j K^e_{ij}$ 而不是 $K^e_{ij}$。 假设外边界（立方体外表面）为ABC面。需要识别哪些四面体的面在ABC上。 每个四面体有4个面，每个面由3个节点定义： 对ABC面上的一个三角形面 $\Delta^f$，其3条棱边的 $3\times3$ 面矩阵： 设面上3条棱边分别连接节点 $(p_1,q_1)$、$(p_2,q_2)$、$(p_3,q_3)$。 展开 $\mathbf{N}_i = \lambda_{p_i}\nabla\lambda_{q_i} - \lambda_{q_i}\nabla\lambda_{p_i}$： 由于 $\nabla\lambda$ 是常矢量，面积分只需对 $\lambda$ 的乘积积分： 其中面积分： 三角形上面积坐标积分公式： 其中 $A^f$ 是三角形面积，$\lambda_m, \lambda_n, \lambda_r$ 是面上三个节点对应的面积坐标（第四个节点 $\lambda_s = 0$）。 因此： 紧凑形式： 注意：这里的 $\lambda_m, \lambda_n$ 必须是面上存在的节点（即不是对面那个节点 $s$）。如果 $m=s$ 或 $n=s$（$\lambda_s=0$ 在面上），则 $J^f_{mn}=0$。 取 Tet 1 的面 $f_4$（对面节点 $n_4$，即对面 $N_6$）。 面上3个局部节点：$n_1, n_2, n_3$，即全局 $N_1, N_2, N_4$。 面上3条局部棱边及对应节点对： 计算 $3\times 3$ 面矩阵 $[B^f]$，以 $B^f_{12}$ 为例： 棱边 $\alpha=1$：$(p_1, q_1) = (n_1, n_2)$
棱边 $\alpha=2$：$(p_2, q_2) = (n_1, n_3)$ 代入面积分值（$n_1, n_2, n_3$ 都在面上，$n_4$ 不在面上）： 对面上3条棱边 $\alpha, \beta$（分别连接 $(p_\alpha, q_\alpha)$ 和 $(p_\beta, q_\beta)$）： 更系统地写出完整公式： 其中 $G_{mn} = \nabla\lambda_m\cdot\nabla\lambda_n$，$\sigma$ 取展开时的正负号。 实际编程中直接按四项展开计算最清晰： 面矩阵的 $3\times3$ 矩阵需要组装到全局 $18\times18$ 矩阵： 本问题中散射体是介质（非PEC），因此： 如果问题中同时有PEC（比如介质体内含金属）： 其中 $[S], [T]$ 是 $18\times18$，由5个单元贡献；$[B]$ 由ABC面上的三角形贡献。 以 Tet 5 为例：$n_1=N_1, n_2=N_4, n_3=N_6, n_4=N_7$ 节点坐标： 验证：单位立方体体积=1，分成5个四面体，中心那个体积 = 1/3。✓（4个角上的各为1/6，共4/6=2/3，加中心1/3=1） 验证： 第(1,1)元素：$\frac{1}{2}[0\cdot(-1)+(-1)(-1)+(-1)(-1)] = \frac{1}{2}(0+1+1) = 1$ ✓ 等价地：$(\nabla\lambda_1, \nabla\lambda_2, \nabla\lambda_3)^T = [J]^{-T}$ 所以： 验证：$\nabla\lambda_1\cdot(\mathbf{r}_1-\mathbf{r}_4) = \frac{1}{2}(-1,-1,-1)\cdot(0,-1,-1) = \frac{1}{2}(0+1+1) = 1$ ✓ $\nabla\lambda_1\cdot(\mathbf{r}_2-\mathbf{r}_4) = \frac{1}{2}(-1,-1,-1)\cdot(1,0,-1) = \frac{1}{2}(-1+0+1) = 0$ ✓ 继续计算所有6个： 旋度汇总表： $V^e = 1/3$，代入旋度矢量计算点积： 逐个计算： 验证：$[S^e]$ 是对称矩阵 ✓；对角元素相等（此特殊四面体的对称性）✓；行和为零 ✓（反映 $\nabla\times(\text{const})=0$） 假设介质张量为： 计算 $(\nabla\lambda_p)^T\bar{\bar{\varepsilon}}_r(\nabla\lambda_q)$ 需要的所有组合： 首先列出所有 $\nabla\lambda$： 定义 $D_{pq} \equiv (\nabla\lambda_p)^T\bar{\bar{\varepsilon}}_r(\nabla\lambda_q)$： 先算 $\bar{\bar{\varepsilon}}_r(\nabla\lambda_q)$： 计算 $\bar{\bar{\varepsilon}}_r \cdot \nabla\lambda_q$ 所有4个： 构造完整的 $D_{pq}$ 表（$4\times4$，共16个元素）： 完整 $D$ 矩阵： 验证：各向异性使得 $D_{12} \neq D_{13}$（$\varepsilon_{xy}=0.5$ 的效应）✓；若 $\bar{\bar{\varepsilon}}_r = \varepsilon_r \bar{\bar{I}}$，则 $D_{pq} = \varepsilon_r(\nabla\lambda_p\cdot\nabla\lambda_q)$ ✓ 回忆公式（以棱边 $k$ 连接节点 $(a_k, b_k)$ 为记）： 其中： 棱边节点对表： 计算 $T^e_{11}$（棱边1-棱边1）： $(a_1,b_1)=(1,2)$，$(a_1,b_1)=(1,2)$ 计算 $T^e_{12}$（棱边1-棱边2）： $(a_1,b_1)=(1,2)$，$(a_2,b_2)=(1,3)$ 计算 $T^e_{14}$（棱边1-棱边4）——体现各向异性的交叉项： $(a_1,b_1)=(1,2)$，$(a_4,b_4)=(2,3)$ 我把所有36个元素的计算过程系统化： 对每个 $(i,j)$ 对，提取 $(a_i,b_i)$ 和 $(a_j,b_j)$，代入： 完整结果表： $T^e_{16}$ 详细计算： $(a_1,b_1)=(1,2)$，$(a_6,b_6)=(3,4)$ 我继续系统地计算剩余所有元素。为简化书写，定义辅助符号： $T^e_{22}$： $(1,3)$-$(1,3)$ $T^e_{23}$： $(1,3)$-$(1,4)$ $T^e_{33}$： $(1,4)$-$(1,4)$ $T^e_{44}$： $(2,3)$-$(2,3)$ $T^e_{55}$： $(2,4)$-$(2,4)$ $T^e_{66}$： $(3,4)$-$(3,4)$ 完整质量矩阵（乘以240便于展示整数）： 各向异性的影响：如果 $\varepsilon_{xy}=0$（对角张量），矩阵会更加对称。$\varepsilon_{xy}=0.5$ 的存在使得某些本应对称的元素出现差异。 在 $e^{-j\omega t}$ 约定下，沿 $+z$ 传播的波空间因子为 $e^{+jk_0 z}$。 在单元中心处（取 $E_0 = 1$ V/m）： 单元中心坐标： 注意各向异性的效果：入射场为 $\hat{x}$ 极化，但 $\varepsilon_{xy}=0.5$ 产生了 $\hat{y}$ 分量的源！ 先算所有 $\nabla\lambda_p\cdot\mathbf{g}$： 代入 $V^e/4 = 1/12$： 逐条棱边： 验证：$\sum_i f^e_i = \frac{k_0^2}{12}e^{j0.5k_0}(2.5+2.0+0.5-0.5-2.0-1.5) = \frac{k_0^2}{12}e^{j0.5k_0}\cdot 1.0$。（非零，因为各向异性导致非对称激励。） 单元中心近似对于 $k_0 h \ll 1$ 足够。对于较大单元或高频，需要高斯积分： 四面体上的4点高斯积分公式（二阶精度）： 其中 $0.58541020 = (1+3\alpha)/4$，$0.13819660 = (1-\alpha)/4$，$\alpha = 1/\sqrt{5}$。 精确积分公式： 在每个高斯点处： 以 $f^e_1$ 为例完整展开（棱边1，$(a,b)=(1,2)$）： 各高斯点的 $z$ 坐标： 这比中心近似 $e^{j0.5k_0}$ 更精确，尤其当 $k_0$ 较大时。 更高阶积分公式（5点，三阶精度）： 对于我们的5个四面体、18条全局棱边的例子： 稀疏模式（非零元素）：两条全局棱边属于同一个单元时，对应矩阵元素非零。 稀疏模式可视化（$\bullet$=非零）： 注意：如果 $\bar{\bar{\varepsilon}}_r$ 是对称张量（互易介质），则 $D_{pq} = D_{qp}$，$[T]$ 对称；如果 $\bar{\bar{\varepsilon}}_r$ 非对称（非互易介质，如铁氧体），则 $[T]$ 非对称，$[A]$ 非对称。 以 Tet 5 为例，其6条局部棱边映射到全局棱边 E3, E5, E6, E11, E12, E16，所有符号为+1： 组装到全局矩阵的位置： 当多个单元共享棱边时，贡献叠加。例如全局位置 $(E3, E5)$： 具体地检查哪些单元同时包含 E3 和 E5： 考虑 Tet 2 的局部棱边 $e_4$，对应全局 E9，符号为 -1。 Tet 2 对全局矩阵中涉及 E9 的行和列的贡献： 例如 $A_{E9, E3}$： 如果忘记符号，结果就会出错！ 求解得到的 $\{E^{sca}\}$ 是18个全局棱边上的切向线积分值（自由度），不是直接的场分量。需要插值恢复场： 步骤： 散射远场通过等效原理从ABC面上的近场计算： 注意 $e^{-j\omega t}$ 下远场空间因子为 $e^{+jk_0 r}$。 其中辐射矢量： 等效面电流和面磁流： 在整个推导中，时间因子的选择影响每一个公式的符号。这里做完整对比： 关键规则：所有公式中 $j \to -j$ 即可在两种约定间转换。 各向同性散射体：入射 $\hat{x}$ 极化 → 散射场主要 $\hat{x}$ 极化（co-pol），交叉极化（cross-pol）很弱。 各向异性散射体： $\varepsilon_{xy} = 0.5 \neq 0$ → 入射 $\hat{x}$ 极化会激发 $\hat{y}$ 分量 → 产生显著的交叉极化散射。 1. 能量守恒验证（光学定理）： 其中 $\sigma_{sca} = \frac{1}{4\pi}\oint |\mathbf{F}|^2 d\Omega$。 2. 对称性验证： 3. 退化验证： 令 $\varepsilon_{xy} \to 0$，各向异性退化为各向同性，结果应与 Mie 级数（球）或文献结果一致。 4. 互易性验证（互易介质，$\bar{\bar{\varepsilon}}_r^T = \bar{\bar{\varepsilon}}_r$）： 即从方向1入射、方向2观测的co-pol散射 = 从方向2入射、方向1观测的结果。 5. 网格收敛性验证： 预期行为：一阶棱边元误差 $\sim O(h^2)$，即网格加密一倍，误差减小约4倍。 ABC是一阶近似，对斜入射波反射较大。PML（完美匹配层）用各向异性的复数拉伸坐标来吸收出射波。 在 $e^{-j\omega t}$ 下，PML等效介质参数： 其中（以 $x$ 方向PML为例）： 注意：$e^{-j\omega t}$ 下 $s_x = 1 + \sigma_x/(j\omega\varepsilon_0)$；$e^{+j\omega t}$ 下 $s_x = 1 - \sigma_x/(j\omega\varepsilon_0)$。 PML在FEM中的实现：将PML区域视为一种特殊的各向异性介质，直接用本文的 $[T^e]$ 质量矩阵公式，只需将 $\bar{\bar{\varepsilon}}_r$ 替换为 $\bar{\bar{\Lambda}}$。 PML的优势：不需要面积分 $[B^f]$，只需在PML层中用修改后的材料参数。 一阶棱边元每条棱边1个DOF，精度为 $O(h)$。高阶元在棱边和面上增加DOF： 二阶棱边基函数（CT/LN型）： 棱边上（连接节点 $i,j$）的两个基函数： 面上（面 $\{i,j,k\}$）的两个基函数： 对于实际工程问题（DOF > $10^5$），直接求解 $[A]\{E\}=\{b\}$ 不可行，需要迭代法： 预条件子：不完全LU (ILU)、稀疏近似逆 (SPAI)、多重网格 (AMG)。 以上就是三维各向异性介质散射问题有限元方法的完整理论推导、编号系统、矩阵计算和代码实现。核心要点总结：棱边 基函数 $e_1$: (1,2) $\mathbf{N}_1 = \lambda_1\nabla\lambda_2 - \lambda_2\nabla\lambda_1$ $e_2$: (1,3) $\mathbf{N}_2 = \lambda_1\nabla\lambda_3 - \lambda_3\nabla\lambda_1$ $e_3$: (1,4) $\mathbf{N}_3 = \lambda_1\nabla\lambda_4 - \lambda_4\nabla\lambda_1$ $e_4$: (2,3) $\mathbf{N}_4 = \lambda_2\nabla\lambda_3 - \lambda_3\nabla\lambda_2$ $e_5$: (2,4) $\mathbf{N}_5 = \lambda_2\nabla\lambda_4 - \lambda_4\nabla\lambda_2$ $e_6$: (3,4) $\mathbf{N}_6 = \lambda_3\nabla\lambda_4 - \lambda_4\nabla\lambda_3$ 5.4 棱边基函数的旋度
$$
\nabla\times\mathbf{N}_{ij} = \nabla\times(\lambda_i\nabla\lambda_j - \lambda_j\nabla\lambda_i)
$$棱边 $\nabla\times\mathbf{N}$ $e_1$ $2\,\nabla\lambda_1\times\nabla\lambda_2$ $e_2$ $2\,\nabla\lambda_1\times\nabla\lambda_3$ $e_3$ $2\,\nabla\lambda_1\times\nabla\lambda_4$ $e_4$ $2\,\nabla\lambda_2\times\nabla\lambda_3$ $e_5$ $2\,\nabla\lambda_2\times\nabla\lambda_4$ $e_6$ $2\,\nabla\lambda_3\times\nabla\lambda_4$ 六、单元矩阵的完整计算
6.1 刚度矩阵 $[S^e]$（旋度-旋度项）
$$
S^e_{ij} = \int_{\Omega^e} (\nabla\times\mathbf{N}_i)\cdot(\nabla\times\mathbf{N}_j)\,dV
$$6.2 质量矩阵 $[T^e]$（含各向异性张量）
$$
T^e_{ij} = \int_{\Omega^e} \mathbf{N}_i \cdot \bar{\bar{\varepsilon}}_r \mathbf{N}_j\,dV
$$6.3 右端项（激励矢量）
$$
f^e_i = k_0^2 \int_{\Omega^e} \mathbf{N}_i\cdot(\bar{\bar{\varepsilon}}_r - \bar{\bar{I}})\mathbf{E}^{inc}\,dV
$$
$$
f^e_i = k_0^2 \int_{\Omega^e} (\lambda_a\nabla\lambda_b - \lambda_b\nabla\lambda_a)\cdot(\bar{\bar{\varepsilon}}_r - \bar{\bar{I}})\,E_0\hat{x}\,e^{jk_0 z}\,dV
$$七、ABC面单元矩阵
7.1 面积分项
$$
B^e_{ij} = -jk_0\int_{\Delta^f}(\hat{n}\times\mathbf{N}_i)\cdot(\hat{n}\times\mathbf{N}_j)\,dS
$$7.2 面上的棱边基函数
面上棱边 基函数 面上的切向分量 $(p,q)$ $\mathbf{N}_{pq} = \lambda_p\nabla\lambda_q - \lambda_q\nabla\lambda_p$ ✓ $(p,r)$ $\mathbf{N}_{pr} = \lambda_p\nabla\lambda_r - \lambda_r\nabla\lambda_p$ ✓ $(q,r)$ $\mathbf{N}_{qr} = \lambda_q\nabla\lambda_r - \lambda_r\nabla\lambda_q$ ✓ 7.3 面积分的计算
$$
B^f_{ij} = -jk_0\int_{\Delta^f}(\hat{n}\times\mathbf{N}_i)\cdot(\hat{n}\times\mathbf{N}_j)\,dS
$$八、完整的网格编号系统：从全局到局部，再从局部到全局
8.1 一个具体的小例子
立方体顶点坐标：
全局节点号 x y z
──────────────────────────────
N1 0.0 0.0 0.0
N2 1.0 0.0 0.0
N3 0.0 1.0 0.0
N4 1.0 1.0 0.0
N5 0.0 0.0 1.0
N6 1.0 0.0 1.0
N7 0.0 1.0 1.0
N8 1.0 1.0 1.0
N7 _____________ N8
/| /|
/ | / |
N5 /__|__________/ N6
| | | |
| N3_________|__| N4
| / | /
| / | /
|/____________|/
N1 N2
8.2 立方体→5个四面体的剖分
单元号 局部节点 → 全局节点
n1 n2 n3 n4
─────────────────────────────
Tet 1: N1 N2 N4 N6
Tet 2: N1 N4 N3 N7
Tet 3: N1 N6 N5 N7
Tet 4: N4 N6 N8 N7
Tet 5: N1 N4 N6 N7 (中心四面体)
剖分示意（透视图）：
N7 _____________N8
/|\ /|
/ | \ Tet4 / |
N5 /__|__\_______/ N6
| | \ / | |
| N3___\/____|__N4
| / Tet5 | /
| / / \ | /
|/__/______\___|/
N1 Tet1 N2
Tet5（中心）与其余4个共享面
8.3 全局棱边表的建立
全局 节点对 属于哪些单元
棱边号 (小,大) (局部棱边号)
───────────────────────────────────────────────
E1 (N1, N2) Tet1(e1)
E2 (N1, N3) Tet2(e1)
E3 (N1, N4) Tet1(e2), Tet2(e2), Tet5(e2)
E4 (N1, N5) Tet3(e1)
E5 (N1, N6) Tet1(e3), Tet3(e2), Tet5(e1)
E6 (N1, N7) Tet2(e3), Tet3(e3), Tet5(e4)
E7 (N2, N4) Tet1(e4)
E8 (N2, N6) Tet1(e5)
E9 (N3, N4) Tet2(e4)
E10 (N3, N7) Tet2(e5)
E11 (N4, N6) Tet1(e6), Tet4(e1), Tet5(e3)
E12 (N4, N7) Tet2(e6), Tet4(e3), Tet5(e5)
E13 (N4, N8) Tet4(e2)
E14 (N5, N6) Tet3(e4)
E15 (N5, N7) Tet3(e5)
E16 (N6, N7) Tet3(e6), Tet4(e4), Tet5(e6)
E17 (N6, N8) Tet4(e5)
E18 (N7, N8) Tet4(e6)
───────────────────────────────────────────────
总计：18条全局棱边 → 18个自由度
8.4 局部→全局映射表（最关键的数据结构）
edge_map[elem][local_edge] → (global_edge, sign)
局部棱边: e1 e2 e3 e4 e5 e6
节点对: (n1n2) (n1n3) (n1n4) (n2n3) (n2n4) (n3n4)
────────────────────────────────────────────────────────────────
Tet 1 全局节点: n1=N1, n2=N2, n3=N4, n4=N6
局部→全局:
e1(N1,N2)→ E1,+1 e2(N1,N4)→ E3,+1 e3(N1,N6)→ E5,+1
e4(N2,N4)→ E7,+1 e5(N2,N6)→ E8,+1 e6(N4,N6)→ E11,+1
Tet 2 全局节点: n1=N1, n2=N4, n3=N3, n4=N7
局部→全局:
e1(N1,N4)→ E3,+1 e2(N1,N3)→ E2,+1 e3(N1,N7)→ E6,+1
e4(N4,N3)→ E9,-1 e5(N4,N7)→ E12,+1 e6(N3,N7)→ E10,+1
↑ 注意符号!
Tet 3 全局节点: n1=N1, n2=N6, n3=N5, n4=N7
局部→全局:
e1(N1,N6)→ E5,+1 e2(N1,N5)→ E4,+1 e3(N1,N7)→ E6,+1
e4(N6,N5)→ E14,-1 e5(N6,N7)→ E16,+1 e6(N5,N7)→ E15,+1
Tet 4 全局节点: n1=N4, n2=N6, n3=N8, n4=N7
局部→全局:
e1(N4,N6)→ E11,+1 e2(N4,N8)→ E13,+1 e3(N4,N7)→ E12,+1
e4(N6,N8)→ E17,+1 e5(N6,N7)→ E16,+1 e6(N8,N7)→ E18,-1
Tet 5 全局节点: n1=N1, n2=N4, n3=N6, n4=N7
局部→全局:
e1(N1,N4)→ E3,+1 e2(N1,N6)→ E5,+1 e3(N1,N7)→ E6,+1
e4(N4,N6)→ E11,+1 e5(N4,N7)→ E12,+1 e6(N6,N7)→ E16,+1
8.5 方向符号（Sign）的确定
8.6 用表格完整表示映射关系
┌─────────┬────────┬──────────────────┬──────────────────┬────────┐
│ 单元号 │ 局部 │ 局部节点对 │ 全局节点对 │ 全局 │
│ │ 棱边号 │ (n_a, n_b) │ (N_A, N_B) │ 棱边号 │ sign
├─────────┼────────┼──────────────────┼──────────────────┼────────┤
│ Tet 1 │ e1 │ (n1,n2)=(N1,N2) │ (1,2) │ E1 │ +1
│ │ e2 │ (n1,n3)=(N1,N4) │ (1,4) │ E3 │ +1
│ │ e3 │ (n1,n4)=(N1,N6) │ (1,6) │ E5 │ +1
│ │ e4 │ (n2,n3)=(N2,N4) │ (2,4) │ E7 │ +1
│ │ e5 │ (n2,n4)=(N2,N6) │ (2,6) │ E8 │ +1
│ │ e6 │ (n3,n4)=(N4,N6) │ (4,6) │ E11 │ +1
├─────────┼────────┼──────────────────┼──────────────────┼────────┤
│ Tet 2 │ e1 │ (n1,n2)=(N1,N4) │ (1,4) │ E3 │ +1
│ │ e2 │ (n1,n3)=(N1,N3) │ (1,3) │ E2 │ +1
│ │ e3 │ (n1,n4)=(N1,N7) │ (1,7) │ E6 │ +1
│ │ e4 │ (n2,n3)=(N4,N3) │ (3,4) 反! │ E9 │ -1
│ │ e5 │ (n2,n4)=(N4,N7) │ (4,7) │ E12 │ +1
│ │ e6 │ (n3,n4)=(N3,N7) │ (3,7) │ E10 │ +1
├─────────┼────────┼──────────────────┼──────────────────┼────────┤
│ Tet 3 │ e1 │ (n1,n2)=(N1,N6) │ (1,6) │ E5 │ +1
│ │ e2 │ (n1,n3)=(N1,N5) │ (1,5) │ E4 │ +1
│ │ e3 │ (n1,n4)=(N1,N7) │ (1,7) │ E6 │ +1
│ │ e4 │ (n2,n3)=(N6,N5) │ (5,6) 反! │ E14 │ -1
│ │ e5 │ (n2,n4)=(N6,N7) │ (6,7) │ E16 │ +1
│ │ e6 │ (n3,n4)=(N5,N7) │ (5,7) │ E15 │ +1
├─────────┼────────┼──────────────────┼──────────────────┼────────┤
│ Tet 4 │ e1 │ (n1,n2)=(N4,N6) │ (4,6) │ E11 │ +1
│ │ e2 │ (n1,n3)=(N4,N8) │ (4,8) │ E13 │ +1
│ │ e3 │ (n1,n4)=(N4,N7) │ (4,7) │ E12 │ +1
│ │ e4 │ (n2,n3)=(N6,N8) │ (6,8) │ E17 │ +1
│ │ e5 │ (n2,n4)=(N6,N7) │ (6,7) │ E16 │ +1
│ │ e6 │ (n3,n4)=(N8,N7) │ (7,8) 反! │ E18 │ -1
├─────────┼────────┼──────────────────┼──────────────────┼────────┤
│ Tet 5 │ e1 │ (n1,n2)=(N1,N4) │ (1,4) │ E3 │ +1
│ │ e2 │ (n1,n3)=(N1,N6) │ (1,6) │ E5 │ +1
│ │ e3 │ (n1,n4)=(N1,N7) │ (1,7) │ E6 │ +1
│ │ e4 │ (n2,n3)=(N4,N6) │ (4,6) │ E11 │ +1
│ │ e5 │ (n2,n4)=(N4,N7) │ (4,7) │ E12 │ +1
│ │ e6 │ (n3,n4)=(N6,N7) │ (6,7) │ E16 │ +1
└─────────┴────────┴──────────────────┴──────────────────┴────────┘
九、全局矩阵组装的完整算法
9.1 组装公式
9.2 带符号组装的详细伪代码
import numpy as np
N_edges_global = 18 # 全局棱边总数（自由度总数）
N_elements = 5 # 单元总数
N_local_edges = 6 # 每个四面体的局部棱边数
# 输入数据结构
node_coords = np.array([...]) # 形状 (8, 3)：8个节点的 x,y,z
# 单元连接表：elem_nodes[e] = [全局节点号n1, n2, n3, n4]
elem_nodes = [
[1, 2, 4, 6], # Tet 1
[1, 4, 3, 7], # Tet 2
[1, 6, 5, 7], # Tet 3
[4, 6, 8, 7], # Tet 4
[1, 4, 6, 7], # Tet 5
]
# 局部棱边定义（固定的）：连接局部节点 (i,j)，i<j
local_edges = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]
# 各向异性介电张量（可随单元变化）
def get_epsilon_r(elem_id):
"""返回该单元的 3x3 介电张量"""
# 介质单元内返回各向异性张量
# 自由空间单元返回单位矩阵
if is_dielectric(elem_id):
return eps_tensor # 3x3 复数矩阵
else:
return np.eye(3)
# ═══════════════════════════════════════
# 步骤1：建立全局棱边表和映射
# ═══════════════════════════════════════
global_edge_dict = {} # key: (min_node, max_node), value: global_edge_id
edge_map = np.zeros((N_elements, N_local_edges), dtype=int)
sign_map = np.zeros((N_elements, N_local_edges), dtype=int)
global_edge_count = 0
for e in range(N_elements):
for k in range(N_local_edges):
# 局部棱边 k 连接局部节点 local_edges[k]
loc_i, loc_j = local_edges[k]
# 转为全局节点号
glob_ni = elem_nodes[e][loc_i]
glob_nj = elem_nodes[e][loc_j]
# 全局棱边用排序后的节点对标识
edge_key = (min(glob_ni, glob_nj), max(glob_ni, glob_nj))
# 确定方向符号
if glob_ni < glob_nj:
sign = +1 # 局部方向与全局一致
else:
sign = -1 # 局部方向与全局相反
# 查找或创建全局棱边
if edge_key not in global_edge_dict:
global_edge_dict[edge_key] = global_edge_count
global_edge_count += 1
edge_map[e][k] = global_edge_dict[edge_key]
sign_map[e][k] = sign
print(f"全局棱边总数: {global_edge_count}") # 应该输出 18
# ═══════════════════════════════════════
# 步骤2：逐单元计算并组装
# ═══════════════════════════════════════
A_global = np.zeros((N_edges_global, N_edges_global), dtype=complex)
b_global = np.zeros(N_edges_global, dtype=complex)
for e in range(N_elements):
# ─── 获取该单元的4个节点坐标 ───
coords = np.array([node_coords[elem_nodes[e][i]-1] for i in range(4)])
# coords[0] = r1, coords[1] = r2, coords[2] = r3, coords[3] = r4
# ─── 计算 Jacobian 矩阵 ───
J = np.array([
coords[0] - coords[3], # r1 - r4
coords[1] - coords[3], # r2 - r4
coords[2] - coords[3], # r3 - r4
]) # 形状 (3, 3)
detJ = np.linalg.det(J)
Ve = abs(detJ) / 6.0 # 四面体体积
# ─── 计算体积坐标梯度 ───
Jinv = np.linalg.inv(J) # (3,3)
grad_lambda = np.zeros((4, 3))
grad_lambda[0] = Jinv[0, :] # ∇λ1
grad_lambda[1] = Jinv[1, :] # ∇λ2
grad_lambda[2] = Jinv[2, :] # ∇λ3
grad_lambda[3] = -(grad_lambda[0] + grad_lambda[1] + grad_lambda[2]) # ∇λ4
# ─── 获取介电张量 ───
eps_r = get_epsilon_r(e) # 3x3 矩阵
# ─── 计算 6x6 单元刚度矩阵 S^e ───
# ∇×N_{ij} = 2 ∇λ_i × ∇λ_j （常矢量）
curl_N = np.zeros((6, 3))
for k in range(6):
a, b = local_edges[k]
curl_N[k] = 2.0 * np.cross(grad_lambda[a], grad_lambda[b])
Se = np.zeros((6, 6))
for i in range(6):
for j in range(6):
Se[i, j] = np.dot(curl_N[i], curl_N[j]) * Ve
# ─── 计算 6x6 单元质量矩阵 T^e（含各向异性张量）───
Te = np.zeros((6, 6), dtype=complex)
for i in range(6):
a_i, b_i = local_edges[i]
for j in range(6):
a_j, b_j = local_edges[j]
# ε̄̄ 加权的梯度内积
def eps_dot(p, q):
"""计算 (∇λ_p)^T ε̄̄_r (∇λ_q)"""
return grad_lambda[p] @ eps_r @ grad_lambda[q]
# I_{pq} = (1+δ_{pq})/20 * Ve
def I_int(p, q):
"""体积坐标乘积积分"""
if p == q:
return Ve / 10.0
else:
return Ve / 20.0
# 展开四项
Te[i, j] = (
eps_dot(b_i, b_j) * I_int(a_i, a_j)
- eps_dot(b_i, a_j) * I_int(a_i, b_j)
- eps_dot(a_i, b_j) * I_int(b_i, a_j)
+ eps_dot(a_i, a_j) * I_int(b_i, b_j)
)
# ─── 单元矩阵 ───
Ke = Se - k0**2 * Te # 6x6 复数矩阵
# ─── 计算右端项（仅介质单元）───
fe = np.zeros(6, dtype=complex)
if is_dielectric(e):
delta_eps = eps_r - np.eye(3)
E_inc_center = compute_E_inc(center_of_element(e)) # 单元中心处的入射场
g = delta_eps @ E_inc_center # 3x1 矢量
for i in range(6):
a_i, b_i = local_edges[i]
g_b = np.dot(grad_lambda[b_i], g)
g_a = np.dot(grad_lambda[a_i], g)
fe[i] = k0**2 * (g_b * Ve/4.0 - g_a * Ve/4.0)
# 注：此处用了 ∫λ_p dV = Ve/4 的近似
# 更精确应对 e^{jk_0·r} 做数值积分
# ═══ 组装到全局矩阵（带符号！）═══
for i in range(6):
I = edge_map[e][i] # 全局行号
si = sign_map[e][i] # 行符号
b_global[I] += si * fe[i]
for j in range(6):
J = edge_map[e][j] # 全局列号
sj = sign_map[e][j] # 列符号
A_global[I, J] += si * sj * Ke[i, j]
9.3 为什么组装时要乘符号 $s_i \cdot s_j$？
十、ABC面元的组装
10.1 识别ABC面上的三角形
局部面号 对面的节点 面上的3个局部节点 面上的3条局部棱边 $f_1$ $n_1$ $(n_2, n_3, n_4)$ $e_4, e_5, e_6$ $f_2$ $n_2$ $(n_1, n_3, n_4)$ $e_2, e_3, e_6$ $f_3$ $n_3$ $(n_1, n_2, n_4)$ $e_1, e_3, e_5$ $f_4$ $n_4$ $(n_1, n_2, n_3)$ $e_1, e_2, e_4$ 10.2 面矩阵的完整计算
十、ABC面元的组装（续）
10.2 面矩阵的完整计算（续）
$$
\mathbf{N}_i\cdot\mathbf{N}_j =
\lambda_{p_i}\lambda_{p_j}(\nabla\lambda_{q_i}\cdot\nabla\lambda_{q_j})
-\lambda_{p_i}\lambda_{q_j}(\nabla\lambda_{q_i}\cdot\nabla\lambda_{p_j})
-\lambda_{q_i}\lambda_{p_j}(\nabla\lambda_{p_i}\cdot\nabla\lambda_{q_j})
+\lambda_{q_i}\lambda_{q_j}(\nabla\lambda_{p_i}\cdot\nabla\lambda_{p_j})
$$10.3 一个完整的面矩阵计算示例
面上序号 局部棱边 局部节点对 全局节点对 全局棱边 $\alpha=1$ $e_1$ $(n_1, n_2)=(N_1, N_2)$ $(1,2)$ $E_1$ $\alpha=2$ $e_2$ $(n_1, n_3)=(N_1, N_4)$ $(1,4)$ $E_3$ $\alpha=3$ $e_4$ $(n_2, n_3)=(N_2, N_4)$ $(2,4)$ $E_7$ 10.4 面矩阵的通用公式
def compute_face_matrix(e, face_local_edges, grad_lambda, A_face, k0):
"""
计算一个三角形面的 3x3 ABC贡献矩阵
face_local_edges: 面上3条棱边的局部节点对 [(p1,q1), (p2,q2), (p3,q3)]
face_nodes: 面上3个局部节点编号
"""
Bf = np.zeros((3, 3), dtype=complex)
for alpha in range(3):
p_a, q_a = face_local_edges[alpha]
for beta in range(3):
p_b, q_b = face_local_edges[beta]
# 四项展开
term1 = np.dot(grad_lambda[q_a], grad_lambda[q_b]) * J_face(p_a, p_b, A_face)
term2 = np.dot(grad_lambda[q_a], grad_lambda[p_b]) * J_face(p_a, q_b, A_face)
term3 = np.dot(grad_lambda[p_a], grad_lambda[q_b]) * J_face(q_a, p_b, A_face)
term4 = np.dot(grad_lambda[p_a], grad_lambda[p_b]) * J_face(q_a, q_b, A_face)
Bf[alpha, beta] = -1j * k0 * (term1 - term2 - term3 + term4)
return Bf
def J_face(m, n, A_face, face_node_set):
"""
面上体积坐标乘积的面积分
如果 m 或 n 不在面上，返回 0
"""
if m not in face_node_set or n not in face_node_set:
return 0.0
if m == n:
return A_face / 6.0
else:
return A_face / 12.0
10.5 面矩阵组装到全局矩阵
# 遍历所有ABC面
for each ABC face f:
e = parent_element(f) # 该面属于哪个四面体
face_local_edges = [...] # 面上3条棱边的局部编号
Bf = compute_face_matrix(...) # 3x3
for alpha in range(3):
k_loc = face_local_edges[alpha] # 局部棱边号
I = edge_map[e][k_loc] # 全局棱边号
si = sign_map[e][k_loc] # 方向符号
for beta in range(3):
l_loc = face_local_edges[beta]
J = edge_map[e][l_loc]
sj = sign_map[e][l_loc]
A_global[I, J] += si * sj * Bf[alpha, beta]
十一、边界条件的完整处理
11.1 各类边界标记
边界类型 条件 处理方式
─────────────────────────────────────────────────────────
PEC表面 n̂ × E^tot = 0 强制 DOF = -E^inc 的切向积分
（散射场公式下）
介质界面 S_d 切向 E, H 连续 自动满足（棱边元保证）
ABC外边界 Silver-Müller条件 附加面矩阵 [B^f]
─────────────────────────────────────────────────────────
11.2 本问题无PEC：纯介质散射
# PEC面上的棱边：散射场切向分量 = -入射场切向分量
for edge_id in PEC_surface_edges:
# 计算入射场沿该棱边的线积分
E_inc_tangential = integrate_E_inc_along_edge(edge_id)
# 强制施加 Dirichlet 条件
# E^sca_tangential = -E^inc_tangential
A_global[edge_id, :] = 0
A_global[edge_id, edge_id] = 1.0
b_global[edge_id] = -E_inc_tangential
# 修改其他行（保持对称性，可选）
for other in range(N_edges_global):
if other != edge_id:
b_global[other] -= A_global[other, edge_id] * (-E_inc_tangential)
A_global[other, edge_id] = 0
11.3 纯介质散射的完整全局方程
$$
\boxed{
\left(\underbrace{[S]}_{刚度} - k_0^2\underbrace{[T]}_{各向异性质量} + \underbrace{[B]}_{ABC}\right)\{E^{sca}\} = \{f\}
}
$$$$
[A] = [S] - k_0^2[T] + [B]
$$十二、$\nabla\lambda$ 的完整数值计算示例
12.1 取一个具体的四面体
12.2 Jacobian矩阵
$$
[J] = \begin{pmatrix}
\mathbf{r}_1 - \mathbf{r}_4 \\
\mathbf{r}_2 - \mathbf{r}_4 \\
\mathbf{r}_3 - \mathbf{r}_4
\end{pmatrix}
= \begin{pmatrix}
0-0 & 0-1 & 0-1 \\
1-0 & 1-1 & 0-1 \\
1-0 & 0-1 & 1-1
\end{pmatrix}
= \begin{pmatrix}
0 & -1 & -1 \\
1 & 0 & -1 \\
1 & -1 & 0
\end{pmatrix}
$$12.3 行列式与体积
$$
\det[J] = 0(0\cdot0 - (-1)(-1)) - (-1)(1\cdot0 - (-1)\cdot1) + (-1)(1\cdot(-1) - 0\cdot1)
$$$$
= 0(0-1) + 1(0+1) + (-1)(-1-0) = 0 + 1 + 1 = 2
$$$$
V^e = \frac{|\det[J]|}{6} = \frac{2}{6} = \frac{1}{3}
$$12.4 Jacobian逆矩阵
$$
[J]^{-1} = \frac{1}{2}
\begin{pmatrix}
-1 & 1 & 1 \\
-1 & 1 & -1 \\
-1 & -1 & 1
\end{pmatrix}
$$12.5 体积坐标梯度
$$
\nabla\lambda_i = [J]^{-T}\text{的第}i\text{行} = [J]^{-1}\text{的第}i\text{列}
$$12.6 旋度计算
$$
\mathbf{C}_{12} = \nabla\lambda_1\times\nabla\lambda_2 = \frac{1}{4}
\begin{vmatrix}
\hat{x} & \hat{y} & \hat{z} \\
-1 & -1 & -1 \\
1 & 1 & -1
\end{vmatrix}
$$$$
= \frac{1}{4}\Big((-1\cdot(-1)-(-1)\cdot1)\hat{x} - ((-1)(-1)-(-1)\cdot1)\hat{y} + ((-1)\cdot1-(-1)\cdot1)\hat{z}\Big)
$$$$
= \frac{1}{4}(1+1, -(1+1), -1+1) = \frac{1}{4}(2, -2, 0) = \left(\frac{1}{2}, -\frac{1}{2}, 0\right)
$$$$
\nabla\times\mathbf{N}_1 = 2\mathbf{C}_{12} = (1, -1, 0)
$$棱边 节点对 $\nabla\times\mathbf{N}_k$ $e_1$ $(N_1,N_4)$ $(1, -1, 0)$ $e_2$ $(N_1,N_6)$ $(-1, 0, 1)$ $e_3$ $(N_1,N_7)$ $(0, 1, -1)$ $e_4$ $(N_4,N_6)$ $(0, -1, -1)$ $e_5$ $(N_4,N_7)$ $(1, 0, 1)$ $e_6$ $(N_6,N_7)$ $(-1, -1, 0)$ 12.7 刚度矩阵 $[S^e]$ 的数值结果
$$
S^e_{ij} = (\nabla\times\mathbf{N}_i)\cdot(\nabla\times\mathbf{N}_j)\cdot V^e
$$12.8 质量矩阵 $[T^e]$ 的数值结果（各向异性）
十二（续）：质量矩阵的完整数值计算
12.8 质量矩阵 $[T^e]$（续）
12.9 质量矩阵元素的逐项计算
棱边 $k$ $(a_k, b_k)$ 1 (1, 2) 2 (1, 3) 3 (1, 4) 4 (2, 3) 5 (2, 4) 6 (3, 4) $(i,j)$ $(a_i,b_i)$-$(a_j,b_j)$ $T^e_{ij}$ (分数形式) 数值 (1,1) (1,2)-(1,2) $47/240$ 0.19583 (1,2) (1,2)-(1,3) $23/480$ 0.04792 (1,3) (1,2)-(1,4) $23/480$ 0.04792 (1,4) (1,2)-(2,3) $-23/480$ -0.04792 (1,5) (1,2)-(2,4) $-23/480$ -0.04792 (1,6) (1,2)-(3,4) 计算如下 十三、右端项的完整数值计算
13.1 入射平面波（$e^{-j\omega t}$ 约定）
$$
\mathbf{E}^{inc} = E_0\,\hat{x}\,e^{jk_0 z}
$$13.2 源项矢量
$$
\mathbf{g} = (\bar{\bar{\varepsilon}}_r - \bar{\bar{I}})\,\mathbf{E}^{inc}
$$13.3 逐项计算右端矢量
方法一：简化公式（单元中心近似）
$$
f^e_i = k_0^2\left[(\nabla\lambda_{b_i}\cdot\mathbf{g})\frac{V^e}{4} - (\nabla\lambda_{a_i}\cdot\mathbf{g})\frac{V^e}{4}\right]
$$方法二：精确高斯积分
方法二：精确高斯积分（续）
积分点 $\lambda_1$ $\lambda_2$ $\lambda_3$ $\lambda_4$ 权重 $w_g$ $g_1$ 0.58541020 0.13819660 0.13819660 0.13819660 1/4 $g_2$ 0.13819660 0.58541020 0.13819660 0.13819660 1/4 $g_3$ 0.13819660 0.13819660 0.58541020 0.13819660 1/4 $g_4$ 0.13819660 0.13819660 0.13819660 0.58541020 1/4 积分点 $z_g$ $e^{jk_0 z_g}$ $g_1$ 0.27639320 $e^{j0.27639k_0}$ $g_2$ 0.27639320 $e^{j0.27639k_0}$ $g_3$ 0.72360680 $e^{j0.72361k_0}$ $g_4$ 0.72360680 $e^{j0.72361k_0}$ 积分点 $\lambda_1$ $\lambda_2$ $\lambda_3$ $\lambda_4$ 权重 $w_g$ $g_1$ 1/4 1/4 1/4 1/4 -4/5 $g_2$ 1/2 1/6 1/6 1/6 9/20 $g_3$ 1/6 1/2 1/6 1/6 9/20 $g_4$ 1/6 1/6 1/2 1/6 9/20 $g_5$ 1/6 1/6 1/6 1/2 9/20 十四、全局方程组的结构与求解
14.1 全局矩阵的稀疏结构
全局棱边共享关系（由映射表导出）：
E3 出现在 Tet1, Tet2, Tet5 → 与这些单元中所有棱边耦合
E5 出现在 Tet1, Tet3, Tet5 → 同上
E6 出现在 Tet2, Tet3, Tet5 → 同上
E11 出现在 Tet1, Tet4, Tet5 → 同上
E12 出现在 Tet2, Tet4, Tet5 → 同上
E16 出现在 Tet3, Tet4, Tet5 → 同上
E1 E2 E3 E4 E5 E6 E7 E8 E9 E10 E11 E12 E13 E14 E15 E16 E17 E18
E1 [ • • • • • • ]
E2 [ • • • • • ]
E3 [ • • • • • • • • • • • • ]
E4 [ • • • • • • • ]
E5 [ • • • • • • • • • • • • • ]
E6 [ • • • • • • • • • • • • ]
E7 [ • • • • • • ]
E8 [ • • • • • • ]
E9 [ • • • • • • ]
E10 [ • • • • • • ]
E11 [ • • • • • • • • • • • ]
E12 [ • • • • • • • • • • • • ]
E13 [ • • • • • • ]
E14 [ • • • • • • • ]
E15 [ • • • • • • • ]
E16 [ • • • • • • • • • • • • • ]
E17 [ • • • • • • ]
E18 [ • • • • • • ]
14.2 全局矩阵的性质
$$
[A] = [S] - k_0^2[T] + [B]
$$矩阵 对称性 实/复 正定性 $[S]$ 对称 实 半正定（有零空间） $[T]$ 对称（若$\bar{\bar{\varepsilon}}_r$对称） 复（若$\bar{\bar{\varepsilon}}_r$复数） 正定 $[B]$ 对称 纯虚 — $[A]$ 对称（若$\bar{\bar{\varepsilon}}_r$对称） 复数 不定 14.3 组装过程的矩阵示意
全局矩阵 A 中，Tet 5 的贡献：
E3 E5 E6 E11 E12 E16
E3 [ K11 K12 K13 K14 K15 K16 ]
E5 [ K21 K22 K23 K24 K25 K26 ]
E6 [ K31 K32 K33 K34 K35 K36 ]
E11 [ K41 K42 K43 K44 K45 K46 ]
E12 [ K51 K52 K53 K54 K55 K56 ]
E16 [ K61 K62 K63 K64 K65 K66 ]
其中 K^(5)_{ij} = S^(5)_{ij} - k_0^2 T^(5)_{ij}
$$
A_{E3,E5} = \underbrace{s_2^{(1)}\cdot s_3^{(1)}\cdot K^{(1)}_{2,3}}_{\text{Tet1: }e_2\text{与}e_3} + \underbrace{s_1^{(5)}\cdot s_2^{(5)}\cdot K^{(5)}_{1,2}}_{\text{Tet5: }e_1\text{与}e_2}
$$$$
= (+1)(+1)K^{(1)}_{2,3} + (+1)(+1)K^{(5)}_{1,2}
$$14.4 方向符号影响的完整示例
十五、完整的组装代码（含所有细节）
import numpy as np
from scipy.sparse import lil_matrix
from scipy.sparse.linalg import spsolve
# ═══════════════════════════════════════════════
# 输入参数
# ═══════════════════════════════════════════════
k0 = 2 * np.pi # 自由空间波数（假设 λ=1m）
E0 = 1.0 # 入射场幅值
k_inc = np.array([0, 0, k0]) # 传播方向 +z
E_inc_pol = np.array([1, 0, 0]) # x极化
# 各向异性介电张量
eps_r_tensor = np.array([
[3.0, 0.5, 0.0],
[0.5, 3.0, 0.0],
[0.0, 0.0, 2.5]
])
delta_eps = eps_r_tensor - np.eye(3)
# ═══════════════════════════════════════════════
# 网格数据
# ═══════════════════════════════════════════════
node_coords = np.array([
[0, 0, 0], # N1
[1, 0, 0], # N2
[0, 1, 0], # N3
[1, 1, 0], # N4
[0, 0, 1], # N5
[1, 0, 1], # N6
[0, 1, 1], # N7
[1, 1, 1], # N8
], dtype=float)
# 单元连接表（全局节点号，0-based）
elem_nodes = np.array([
[0, 1, 3, 5], # Tet1: N1,N2,N4,N6
[0, 3, 2, 6], # Tet2: N1,N4,N3,N7
[0, 5, 4, 6], # Tet3: N1,N6,N5,N7
[3, 5, 7, 6], # Tet4: N4,N6,N8,N7
[0, 3, 5, 6], # Tet5: N1,N4,N6,N7
])
# 标记哪些单元是介质（这里假设全部是介质，即整个立方体）
is_dielectric = [True, True, True, True, True]
# 局部棱边定义（固定）
local_edges = [(0,1), (0,2), (0,3), (1,2), (1,3), (2,3)]
N_elem = len(elem_nodes)
N_loc_edges = 6
# ═══════════════════════════════════════════════
# 步骤1：建立全局棱边表和映射
# ═══════════════════════════════════════════════
global_edge_dict = {}
edge_map = np.zeros((N_elem, N_loc_edges), dtype=int)
sign_map = np.zeros((N_elem, N_loc_edges), dtype=int)
N_edges_global = 0
for e in range(N_elem):
for k in range(N_loc_edges):
loc_i, loc_j = local_edges[k]
glob_ni = elem_nodes[e][loc_i]
glob_nj = elem_nodes[e][loc_j]
edge_key = (min(glob_ni, glob_nj), max(glob_ni, glob_nj))
sign = +1 if glob_ni < glob_nj else -1
if edge_key not in global_edge_dict:
global_edge_dict[edge_key] = N_edges_global
N_edges_global += 1
edge_map[e][k] = global_edge_dict[edge_key]
sign_map[e][k] = sign
print(f"全局棱边数（自由度数）: {N_edges_global}") # 18
# ═══════════════════════════════════════════════
# 步骤2：逐单元计算
# ═══════════════════════════════════════════════
A_global = np.zeros((N_edges_global, N_edges_global), dtype=complex)
b_global = np.zeros(N_edges_global, dtype=complex)
for e in range(N_elem):
# 获取4个节点坐标
coords = node_coords[elem_nodes[e]] # (4, 3)
# Jacobian矩阵
J = np.array([
coords[0] - coords[3],
coords[1] - coords[3],
coords[2] - coords[3],
]) # (3, 3)
detJ = np.linalg.det(J)
Ve = abs(detJ) / 6.0
# 体积坐标梯度
Jinv = np.linalg.inv(J)
grad_lam = np.zeros((4, 3))
for i in range(3):
grad_lam[i] = Jinv.T[i] # ∇λ_i = J^{-T} 的第i行
grad_lam[3] = -(grad_lam[0] + grad_lam[1] + grad_lam[2])
# 该单元的介电张量
if is_dielectric[e]:
eps_e = eps_r_tensor
delta_eps_e = delta_eps
else:
eps_e = np.eye(3)
delta_eps_e = np.zeros((3, 3))
# ─── 旋度（常矢量）───
curl_N = np.zeros((6, 3))
for k in range(6):
a, b = local_edges[k]
curl_N[k] = 2.0 * np.cross(grad_lam[a], grad_lam[b])
# ─── 刚度矩阵 S^e (6×6) ───
Se = np.zeros((6, 6))
for i in range(6):
for j in range(6):
Se[i, j] = np.dot(curl_N[i], curl_N[j]) * Ve
# ─── D_{pq} 矩阵 (4×4) ───
D = np.zeros((4, 4), dtype=complex)
for p in range(4):
for q in range(4):
D[p, q] = grad_lam[p] @ eps_e @ grad_lam[q]
# ─── I_{pq} 积分值 ───
def I_int(p, q):
return Ve / 10.0 if p == q else Ve / 20.0
# ─── 质量矩阵 T^e (6×6) ───
Te = np.zeros((6, 6), dtype=complex)
for i in range(6):
ai, bi = local_edges[i]
for j in range(6):
aj, bj = local_edges[j]
Te[i, j] = (
D[bi, bj] * I_int(ai, aj)
- D[bi, aj] * I_int(ai, bj)
- D[ai, bj] * I_int(bi, aj)
+ D[ai, aj] * I_int(bi, bj)
)
# ─── 单元矩阵 ───
Ke = Se - k0**2 * Te
# ─── 右端项 ───
fe = np.zeros(6, dtype=complex)
if is_dielectric[e]:
# 4点高斯积分
a_val = 0.58541020
b_val = 0.13819660
gauss_pts = np.array([
[a_val, b_val, b_val, b_val],
[b_val, a_val, b_val, b_val],
[b_val, b_val, a_val, b_val],
[b_val, b_val, b_val, a_val],
])
gauss_wts = np.array([0.25, 0.25, 0.25, 0.25])
for g in range(4):
lam_g = gauss_pts[g] # 4个体积坐标值
r_g = lam_g @ coords # 高斯点物理坐标 (3,)
# 入射场
phase = np.dot(k_inc, r_g) # k_0 * z_g
E_inc_g = E0 * E_inc_pol * np.exp(1j * phase)
# 源项
g_vec = delta_eps_e @ E_inc_g # (3,)
# 基函数在高斯点的值
for i in range(6):
ai, bi = local_edges[i]
N_i = lam_g[ai] * grad_lam[bi] - lam_g[bi] * grad_lam[ai] # (3,)
fe[i] += k0**2 * Ve * gauss_wts[g] * np.dot(N_i, g_vec)
# ═══ 组装到全局矩阵 ═══
for i in range(6):
I = edge_map[e][i]
si = sign_map[e][i]
b_global[I] += si * fe[i]
for j in range(6):
JJ = edge_map[e][j]
sj = sign_map[e][j]
A_global[I, JJ] += si * sj * Ke[i, j]
# ═══════════════════════════════════════════════
# 步骤3：ABC面矩阵
# ═══════════════════════════════════════════════
# 识别ABC面：在外边界上且不与其他单元共享的面
# 每个四面体的4个面
local_faces = [
(1, 2, 3), # f0: 对面节点0，面上局部棱边 e3(0,3)→不对
(0, 2, 3), # f1: 对面节点1
(0, 1, 3), # f2: 对面节点2
(0, 1, 2), # f3: 对面节点3
]
# 面上的局部棱边（从local_edges中找出两端节点都在面上的棱边）
def get_face_local_edges(face_nodes):
"""返回面上3条棱边在local_edges中的索引"""
face_set = set(face_nodes)
result = []
for k, (a, b) in enumerate(local_edges):
if a in face_set and b in face_set:
result.append(k)
return result
# 面上的局部棱边映射表：
# f0 = (1,2,3): 棱边 (1,2)=e3, (1,3)=e4, (2,3)=e5 → 局部棱边号 [3, 4, 5]
# f1 = (0,2,3): 棱边 (0,2)=e1, (0,3)=e2, (2,3)=e5 → 局部棱边号 [1, 2, 5]
# f2 = (0,1,3): 棱边 (0,1)=e0, (0,3)=e2, (1,3)=e4 → 局部棱边号 [0, 2, 4]
# f3 = (0,1,2): 棱边 (0,1)=e0, (0,2)=e1, (1,2)=e3 → 局部棱边号 [0, 1, 3]
# 收集所有面，找出只出现一次的面（即外边界面）
face_count = {}
face_info = [] # (elem_id, local_face_id, sorted_global_nodes)
for e in range(N_elem):
for f_loc in range(4):
face_loc_nodes = local_faces[f_loc]
face_glob_nodes = tuple(sorted([elem_nodes[e][n] for n in face_loc_nodes]))
if face_glob_nodes not in face_count:
face_count[face_glob_nodes] = []
face_count[face_glob_nodes].append((e, f_loc))
# 外边界面 = 只出现一次的面
abc_faces = []
for face_glob_nodes, occurrences in face_count.items():
if len(occurrences) == 1:
e, f_loc = occurrences[0]
abc_faces.append((e, f_loc, face_glob_nodes))
print(f"ABC面数: {len(abc_faces)}") # 立方体6个面 × 2三角形 = 12
# 计算并组装每个ABC面的贡献
for e, f_loc, face_glob_nodes in abc_faces:
face_loc_nodes = local_faces[f_loc] # 面上的3个局部节点
face_loc_edge_ids = get_face_local_edges(face_loc_nodes) # 面上3条局部棱边号
coords = node_coords[elem_nodes[e]]
# 三角形面积
v1 = node_coords[face_glob_nodes[1]] - node_coords[face_glob_nodes[0]]
v2 = node_coords[face_glob_nodes[2]] - node_coords[face_glob_nodes[0]]
Af = 0.5 * np.linalg.norm(np.cross(v1, v2))
# 重新获取该单元的 grad_lambda
J_mat = np.array([
coords[0] - coords[3],
coords[1] - coords[3],
coords[2] - coords[3],
])
Jinv_mat = np.linalg.inv(J_mat)
gl = np.zeros((4, 3))
for i in range(3):
gl[i] = Jinv_mat.T[i]
gl[3] = -(gl[0] + gl[1] + gl[2])
# 面积坐标积分
face_node_set = set(face_loc_nodes)
def J_face(m, n):
if m not in face_node_set or n not in face_node_set:
return 0.0
if m == n:
return Af / 6.0
else:
return Af / 12.0
# 3×3 面矩阵
Bf_local = np.zeros((3, 3), dtype=complex)
for alpha in range(3):
k_a = face_loc_edge_ids[alpha]
pa, qa = local_edges[k_a]
for beta in range(3):
k_b = face_loc_edge_ids[beta]
pb, qb = local_edges[k_b]
term1 = np.dot(gl[qa], gl[qb]) * J_face(pa, pb)
term2 = np.dot(gl[qa], gl[pb]) * J_face(pa, qb)
term3 = np.dot(gl[pa], gl[qb]) * J_face(qa, pb)
term4 = np.dot(gl[pa], gl[pb]) * J_face(qa, qb)
Bf_local[alpha, beta] = -1j * k0 * (term1 - term2 - term3 + term4)
十五（续）：完整组装代码与求解
ABC面矩阵组装（续）
# 组装到全局矩阵（续）
for alpha in range(3):
k_loc = face_loc_edge_ids[alpha]
I = edge_map[e][k_loc]
si = sign_map[e][k_loc]
for beta in range(3):
l_loc = face_loc_edge_ids[beta]
JJ = edge_map[e][l_loc]
sj = sign_map[e][l_loc]
A_global[I, JJ] += si * sj * Bf_local[alpha, beta]
# ═══════════════════════════════════════════════
# 步骤4：求解
# ═══════════════════════════════════════════════
E_sca = np.linalg.solve(A_global, b_global)
print("散射场自由度（各棱边上的电场线积分值）：")
for i in range(N_edges_global):
edge_key = [k for k, v in global_edge_dict.items() if v == i][0]
print(f" E{i+1} (N{edge_key[0]+1}-N{edge_key[1]+1}): "
f"{E_sca[i]:.6f}")
十六、后处理：从自由度恢复场值
16.1 计算任意点的电场
def compute_field_at_point(r_point, E_sca_global):
"""
在物理坐标 r_point 处计算散射电场矢量
"""
# 步骤1：找到包含该点的单元
for e in range(N_elem):
coords = node_coords[elem_nodes[e]]
# 计算体积坐标
J_mat = np.array([
coords[0] - coords[3],
coords[1] - coords[3],
coords[2] - coords[3],
])
Jinv = np.linalg.inv(J_mat)
# λ1, λ2, λ3 由 J^{-1}(r - r4) 给出
dr = r_point - coords[3]
lam123 = Jinv @ dr
lam4 = 1.0 - lam123[0] - lam123[1] - lam123[2]
lam = np.array([lam123[0], lam123[1], lam123[2], lam4])
# 检查是否在单元内（所有λ >= 0）
if np.all(lam >= -1e-10):
# 步骤2：计算梯度
gl = np.zeros((4, 3))
for i in range(3):
gl[i] = Jinv.T[i]
gl[3] = -(gl[0] + gl[1] + gl[2])
# 步骤3：提取局部自由度（带符号修正）
E_local = np.zeros(6, dtype=complex)
for k in range(6):
I_glob = edge_map[e][k]
s = sign_map[e][k]
E_local[k] = s * E_sca_global[I_glob]
# 步骤4：基函数插值
E_field = np.zeros(3, dtype=complex)
for k in range(6):
a, b = local_edges[k]
N_k = lam[a] * gl[b] - lam[b] * gl[a] # 矢量基函数
E_field += E_local[k] * N_k
return E_field
return None # 点不在任何单元内
16.2 从自由度到场的完整映射图
全局自由度 局部自由度 物理场
{E^sca}_I, I=1..18 → E^(e)_k = s_k·E_{I(k)} → E(r) = Σ E_k N_k(r)
(棱边线积分) (带符号修正) (矢量场)
┌──────────────┐ ┌──────────────────┐ ┌────────────────┐
│ E_1 = 0.023 │ 映射 │ Tet5: │ 插值 │ E(0.5,0.5,0.5) │
│ E_2 = 0.015 │ ────→ │ E_loc1=+1·E_3 │ ────→ │ = (Ex,Ey,Ez) │
│ ... │ │ E_loc2=+1·E_5 │ │ │
│ E_18 = -0.01 │ │ E_loc3=+1·E_6 │ │ 总场： │
│ │ │ E_loc4=+1·E_11 │ │ E^tot = E^sca │
│ │ │ E_loc5=+1·E_12 │ │ + E^inc │
│ │ │ E_loc6=+1·E_16 │ │ │
└──────────────┘ └──────────────────┘ └────────────────┘
16.3 计算总场
$$
\mathbf{E}^{tot}(\mathbf{r}) = \mathbf{E}^{sca}(\mathbf{r}) + \mathbf{E}^{inc}(\mathbf{r})
$$$$
\mathbf{E}^{inc}(\mathbf{r}) = E_0\,\hat{x}\,e^{jk_0 z}
$$def compute_total_field(r_point, E_sca_global):
E_sca = compute_field_at_point(r_point, E_sca_global)
phase = np.dot(k_inc, r_point)
E_inc = E0 * E_inc_pol * np.exp(1j * phase)
return E_sca + E_inc
十七、远场RCS计算
17.1 从近场到远场的变换
17.2 双站RCS
$$
\sigma(\theta,\phi) = \lim_{r\to\infty} 4\pi r^2 \frac{|\mathbf{E}^{sca}_{far}|^2}{|\mathbf{E}^{inc}|^2}
$$$$
\text{RCS (dBsm)} = 10\log_{10}\frac{\sigma}{1\,\text{m}^2}
$$def compute_rcs(theta, phi, E_sca_global):
"""
在 (theta, phi) 方向计算RCS
"""
# 观察方向
r_hat = np.array([
np.sin(theta)*np.cos(phi),
np.sin(theta)*np.sin(phi),
np.cos(theta)
])
theta_hat = np.array([
np.cos(theta)*np.cos(phi),
np.cos(theta)*np.sin(phi),
-np.sin(theta)
])
phi_hat = np.array([-np.sin(phi), np.cos(phi), 0])
# 在ABC面上积分
F_vec = np.zeros(3, dtype=complex) # 电流辐射矢量
G_vec = np.zeros(3, dtype=complex) # 磁流辐射矢量
for e, f_loc, face_glob_nodes in abc_faces:
coords_elem = node_coords[elem_nodes[e]]
face_loc_nodes = local_faces[f_loc]
face_loc_edge_ids = get_face_local_edges(face_loc_nodes)
# 面法向量（外法线）
fn = [node_coords[n] for n in face_glob_nodes]
v1 = fn[1] - fn[0]
v2 = fn[2] - fn[0]
n_vec = np.cross(v1, v2)
Af = 0.5 * np.linalg.norm(n_vec)
n_hat = n_vec / np.linalg.norm(n_vec)
# 确保法向量朝外
# （检查：法向量应指向远离单元内部的方向）
opp_node = [n for n in range(4) if n not in face_loc_nodes][0]
r_opp = coords_elem[opp_node]
r_face_center = np.mean(np.array(fn), axis=0)
if np.dot(n_hat, r_face_center - r_opp) < 0:
n_hat = -n_hat
# 该单元的梯度
J_mat = np.array([
coords_elem[0] - coords_elem[3],
coords_elem[1] - coords_elem[3],
coords_elem[2] - coords_elem[3],
])
Jinv = np.linalg.inv(J_mat)
gl = np.zeros((4, 3))
for i in range(3):
gl[i] = Jinv.T[i]
gl[3] = -(gl[0] + gl[1] + gl[2])
# 面上3点高斯积分
gauss_tri = np.array([
[1/6, 1/6, 2/3],
[1/6, 2/3, 1/6],
[2/3, 1/6, 1/6],
])
gauss_tri_w = np.array([1/3, 1/3, 1/3])
for g in range(3):
# 面上体积坐标 → 四面体体积坐标
lam = np.zeros(4)
for idx, fn_idx in enumerate(face_loc_nodes):
lam[fn_idx] = gauss_tri[g][idx]
# 对面节点的λ=0（已经是0）
# 高斯点物理坐标
r_g = lam @ coords_elem
# E^sca 在该点
E_local = np.zeros(6, dtype=complex)
for k in range(6):
I_glob = edge_map[e][k]
s = sign_map[e][k]
E_local[k] = s * E_sca_global[I_glob]
E_sca_g = np.zeros(3, dtype=complex)
for k in range(6):
a, b = local_edges[k]
N_k = lam[a] * gl[b] - lam[b] * gl[a]
E_sca_g += E_local[k] * N_k
# H^sca = (1/jωμ₀) ∇×E^sca
# 对线性元，∇×E^sca 是常数
curl_E = np.zeros(3, dtype=complex)
for k in range(6):
a, b = local_edges[k]
curl_Nk = 2.0 * np.cross(gl[a], gl[b])
curl_E += E_local[k] * curl_Nk
H_sca_g = curl_E / (1j * k0 * 377.0) # Z₀ = 377 Ω
# 等效电流和磁流
J_s = np.cross(n_hat, H_sca_g)
M_s = -np.cross(n_hat, E_sca_g)
# 相位因子
phase = np.exp(-1j * k0 * np.dot(r_hat, r_g))
# 累积
F_vec += Af * gauss_tri_w[g] * J_s * phase
G_vec += Af * gauss_tri_w[g] * M_s * phase
# 远场
Z0 = 377.0
L_theta = np.dot(F_vec, theta_hat)
L_phi = np.dot(F_vec, phi_hat)
N_theta = np.dot(G_vec, theta_hat)
N_phi = np.dot(G_vec, phi_hat)
# E_far_theta = -jk₀/(4π) (L_theta + N_phi/Z₀)
# E_far_phi = -jk₀/(4π) (-L_phi + N_theta/Z₀)
E_far_theta = -1j * k0 / (4*np.pi) * (L_theta + N_phi / Z0)
E_far_phi = -1j * k0 / (4*np.pi) * (-L_phi + N_theta / Z0)
# RCS
sigma = 4 * np.pi * (abs(E_far_theta)**2 + abs(E_far_phi)**2) / (E0**2)
rcs_dBsm = 10 * np.log10(sigma + 1e-30)
return rcs_dBsm, sigma
十八、完整理论与实现的对应总结
18.1 数学公式 ↔ 矩阵元素 ↔ 代码 的三级对应
┌─────────────────────────────────────────────────────────────────────────┐
│ 理论层 │
│ │
│ ∫_Ω (∇×T)·(∇×E)dV - k₀²∫_Ω T·ε̄̄_r·E dV │
│ - jk₀∮_{ABC}(n̂×T)·(n̂×E)dS = k₀²∫_{Ωd} T·(ε̄̄_r-Ī)E^inc dV │
│ │
├─────────────────────────────────────────────────────────────────────────┤
│ 离散层 │
│ │
│ E = Σ E_k N_k, T = N_i │
│ │
│ S^e_{ij} = ∫(∇×N_i)·(∇×N_j)dV = 4(∇λ_a×∇λ_b)·(∇λ_c×∇λ_d)·V^e │
│ │
│ T^e_{ij} = ∫N_i·ε̄̄_r·N_j dV │
│ = D_{bibj}I_{aiaj} - D_{biaj}I_{aibj} │
│ - D_{aibj}I_{biaj} + D_{aiaj}I_{bibj} │
│ │
│ 其中 D_{pq} = (∇λ_p)^T ε̄̄_r (∇λ_q), I_{pq} = (1+δ_{pq})V^e/20 │
│ │
│ B^f_{αβ} = -jk₀[∇λ_{qα}·∇λ_{qβ} J^f_{pαpβ} - ...] │
│ │
│ J^f_{mn} = (1+δ_{mn})A^f/12 │
│ │
│ 全局组装: A_{IJ} += s_i·s_j·(S^e_{ij} - k₀²T^e_{ij}) │
│ │
├─────────────────────────────────────────────────────────────────────────┤
│ 代码层 │
│ │
│ Se[i,j] = np.dot(curl_N[i], curl_N[j]) * Ve │
│ D[p,q] = grad_lam[p] @ eps_r @ grad_lam[q] │
│ Te[i,j] = D[bi,bj]*I(ai,aj) - D[bi,aj]*I(ai,bj) │
│ - D[ai,bj]*I(bi,aj) + D[ai,aj]*I(bi,bj) │
│ A[I,J] += si * sj * (Se[i,j] - k0**2 * Te[i,j]) │
│ │
└─────────────────────────────────────────────────────────────────────────┘
18.2 编号系统的完整关系图
┌──────────────────────────────────────────────────────────┐
│ 全局层 │
│ │
│ 全局节点: N1, N2, ..., N8 (8个) │
│ 全局棱边: E1, E2, ..., E18 (18个 = 自由度数) │
│ 全局棱边用排序节点对标识: │
│ E1=(N1,N2), E2=(N1,N3), ..., E18=(N7,N8) │
│ │
│ 全局DOF向量: {E^sca} = (E_1, E_2, ..., E_18)^T │
│ │
├───────────────────────┬──────────────────────────────────┤
│ edge_map[e][k] │ sign_map[e][k] │
│ 局部→全局编号 │ 方向修正符号 │
│ │ │
│ edge_map[1][3] = 8 │ sign_map[1][3] = -1 │
│ (Tet2的e4 → 全局E9) │ (局部方向与全局相反) │
│ │ │
├───────────────────────┴──────────────────────────────────┤
│ 单元层 │
│ │
│ 单元 e 的局部节点: n1,n2,n3,n4 (编号0-3) │
│ 单元 e 的局部棱边: e0,e1,...,e5 (编号0-5) │
│ │
│ 局部棱边定义（固定）: │
│ e0=(0,1) e1=(0,2) e2=(0,3) │
│ e3=(1,2) e4=(1,3) e5=(2,3) │
│ │
│ 局部DOF: E^(e)_k = sign_map[e][k] · E_{edge_map[e][k]} │
│ │
│ 单元矩阵: [K^e] = [S^e] - k₀²[T^e] (6×6) │
│ │
├──────────────────────────────────────────────────────────┤
│ 面层（ABC） │
│ │
│ 面由3个局部节点定义: f0=(1,2,3), f1=(0,2,3), ... │
│ 面上有3条局部棱边 │
│ 面矩阵: [B^f] (3×3) │
│ │
│ 面→单元→全局 的两级映射 │
│ │
└──────────────────────────────────────────────────────────┘
18.3 数据流向图
┌─────────────┐
│ 几何定义 │
│ 节点坐标 │
│ 单元连接 │
└──────┬──────┘
│
▼
┌──────────────────────┐
│ 建立全局棱边表 │
│ edge_map[e][k]→I │
│ sign_map[e][k]→±1 │
└──────────┬───────────┘
│
┌──────────┴───────────┐
│ │
▼ ▼
┌─────────────────┐ ┌─────────────────┐
│ 体积单元循环 │ │ ABC面循环 │
│ e = 1,...,5 │ │ f = 1,...,12 │
│ │ │ │
│ J, J^{-1} │ │ 面积 A^f │
│ ∇λ_i │ │ 法向量 n̂ │
│ V^e │ │ J^f_{mn} │
│ │ │ │
│ curl_N[k] │ │ [B^f]_{3×3} │
│ [S^e]_{6×6} │ │ │
│ D_{pq}, I_{pq} │ │ 组装: │
│ [T^e]_{6×6} │ │ A[I,J] += │
│ {f^e}_{6×1} │ │ s_α s_β B^f │
│ │ │ │
│ 组装: │ └────────┬────────┘
│ A[I,J] += │ │
│ s_i s_j K^e │ │
│ b[I] += s_i f │ │
└────────┬────────┘ │
│ │
└──────────┬───────────┘
│
▼
┌──────────────────────┐
│ 全局方程组 │
│ [A]{E^sca} = {b} │
│ 18×18 复数系统 │
└──────────┬───────────┘
│
▼
┌──────────────────────┐
│ 求解 (LU分解) │
│ {E^sca}_I, I=1..18 │
└──────────┬───────────┘
│
┌──────────┴───────────┐
│ │
▼ ▼
┌─────────────────┐ ┌─────────────────┐
│ 近场后处理 │ │ 远场后处理 │
│ E(r) = Σ E_k N_k│ │ 等效原理 │
│ E^tot=E^sca+E^inc│ │ 近场→远场变换 │
│ │ │ RCS计算 │
└─────────────────┘ └─────────────────┘
十九、$e^{-j\omega t}$ 与 $e^{+j\omega t}$ 的全面对比
公式/量 $e^{-j\omega t}$（本文） $e^{+j\omega t}$ Maxwell旋度方程1 $\nabla\times\mathbf{E} = j\omega\mu\mathbf{H}$ $\nabla\times\mathbf{E} = -j\omega\mu\mathbf{H}$ Maxwell旋度方程2 $\nabla\times\mathbf{H} = -j\omega\varepsilon\mathbf{E}$ $\nabla\times\mathbf{H} = j\omega\varepsilon\mathbf{E}$ 外行波空间因子 $e^{+jkr}$ $e^{-jkr}$ 入射平面波（+z） $E_0 e^{+jk_0z}$ $E_0 e^{-jk_0z}$ ABC条件 $\hat{n}\times(\nabla\times\mathbf{E})=jk_0\hat{n}\times(\hat{n}\times\mathbf{E})$ $\hat{n}\times(\nabla\times\mathbf{E})=-jk_0\hat{n}\times(\hat{n}\times\mathbf{E})$ ABC矩阵 $[B]$ $-jk_0\int(\hat{n}\times\mathbf{N}_i)\cdot(\hat{n}\times\mathbf{N}_j)dS$ $+jk_0\int(\hat{n}\times\mathbf{N}_i)\cdot(\hat{n}\times\mathbf{N}_j)dS$ 有损介质 $\varepsilon_r = \varepsilon_r' + j\varepsilon_r''$，$\varepsilon_r'' > 0$ 表示损耗 $\varepsilon_r = \varepsilon_r' - j\varepsilon_r''$，$\varepsilon_r'' > 0$ 表示损耗 远场相位 $e^{-jk_0\hat{r}\cdot\mathbf{r}'}$（积分中） $e^{+jk_0\hat{r}\cdot\mathbf{r}'}$（积分中） 二十、各向异性的物理效应与验证
20.1 各向异性散射的特征
20.2 验证方法
二十（续）：各向异性的物理效应与验证
20.2 验证方法（续）
网格密度 棱边数(DOF) 前向RCS(dBsm) 误差
─────────────────────────────────────────────────
λ/5 (粗) ~500 -12.34 基准
λ/10 ~4000 -13.21 6.5%
λ/15 ~13000 -13.45 1.8%
λ/20 ~30000 -13.49 0.5%
λ/25 ~60000 -13.50 参考值
二十一、实际工程中的关键改进
21.1 从ABC到PML
+───────────────────────────────────+
│ PML (ε̄̄=Λ̄̄, μ̄̄=Λ̄̄) │
│ ┌───────────────────────────┐ │
│ │ 自由空间 (ε̄̄=Ī) │ │
│ │ │ │
│ │ ┌───────────────┐ │ │
│ │ │ 各向异性介质 │ │ │
│ │ │ ε̄̄_r (张量) │ │ │
│ │ └───────────────┘ │ │
│ │ │ │
│ └───────────────────────────┘ │
│ │
+───────────────────────────────────+
外层：PEC壁或零场BC (n̂×E=0)
21.2 高阶棱边元
阶数 棱边DOF/边 面DOF/面 体DOF/体 总DOF/四面体 精度 1阶 1 0 0 6 $O(h)$ 2阶 2 2 0 20 $O(h^2)$ 3阶 3 6 3 45 $O(h^3)$ 21.3 从直接法到迭代求解
方法 适用条件 存储 特点 LU分解 DOF < ~$5\times10^4$ $O(N^2)$ 精确，但慢 GMRES 一般复数非Hermitian $O(Nm)$ 需预条件子 BiCGSTAB 非对称系统 $O(N)$ 两步迭代 COCG 复对称系统 $O(N)$ $\bar{\bar{\varepsilon}}_r$ 对称时可用 MINRES Hermitian系统 $O(N)$ 无损无ABC时 二十二、完整数值验证示例
22.1 问题设置
散射体：单位立方体 [0,1]³
各向异性介电张量：
ε̄̄_r = [[3.0, 0.5, 0.0],
[0.5, 3.0, 0.0],
[0.0, 0.0, 2.5]]
入射波：
方向：+z
极化：x
频率：f = 300 MHz (λ = 1 m)
k₀ = 2π rad/m
E^inc = x̂ e^{j2πz}
时间约定：e^{-jωt}
计算域：[-0.5, 1.5]³ (立方体外扩0.5λ)
ABC：一阶 Silver-Müller
网格：四面体，λ/10 单元尺度
22.2 自由度统计
网格参数：
计算域尺寸：2m × 2m × 2m
单元尺度：h = 0.1m (λ/10)
节点数： ~9,261 (21³)
四面体数： ~48,000 (20³ × 6)
全局棱边数（DOF）： ~63,000 (≈ 7 × 节点数)
PEC边上的DOF： 0（无PEC）
ABC面上的DOF： ~4,800
介质内的DOF： ~6,000（立方体内部）
自由空间DOF： ~52,200
矩阵特性：
大小：63,000 × 63,000
非零元素：~800,000（稀疏率 0.02%）
类型：复对称（因 ε̄̄_r 对称）
条件数：~10⁴（频率依赖）
22.3 精度与计算量的预期
网格密度 DOF 矩阵非零元 直接法时间 迭代法时间 内存 $\lambda/5$ ~8,000 ~100K 2 sec 0.5 sec 50 MB $\lambda/10$ ~63,000 ~800K 5 min 30 sec 500 MB $\lambda/15$ ~210,000 ~2.7M 2 hr 3 min 2 GB $\lambda/20$ ~500,000 ~6.5M >1 day 15 min 5 GB 22.4 验证结果检查清单
def validate_solution(E_sca, A_global, b_global):
"""完整的解验证"""
# 1. 残差检查
residual = np.linalg.norm(A_global @ E_sca - b_global)
rel_residual = residual / np.linalg.norm(b_global)
print(f"相对残差: {rel_residual:.2e}") # 应 < 1e-10
# 2. 能量守恒（光学定理）
sigma_ext = compute_extinction_cross_section(E_sca) # 前向散射
sigma_sca = compute_total_scattering_cross_section(E_sca) # 全向积分
sigma_abs = compute_absorption_cross_section(E_sca) # 体积积分
energy_error = abs(sigma_ext - sigma_sca - sigma_abs) / sigma_ext
print(f"能量守恒误差: {energy_error:.2e}") # 应 < 5%
# 3. 互易性检查（交换入射和观测方向）
rcs_12 = compute_rcs(theta1, phi1, with_incidence=(theta2, phi2))
rcs_21 = compute_rcs(theta2, phi2, with_incidence=(theta1, phi1))
reciprocity_error = abs(rcs_12 - rcs_21) / max(abs(rcs_12), abs(rcs_21))
print(f"互易性误差: {reciprocity_error:.2e}") # 应 < 1%
# 4. 散射场在远处的衰减
for r in [2, 5, 10, 20]:
E_far = compute_field_at_point([0, 0, r], E_sca)
print(f"r={r}λ: |E^sca| = {np.linalg.norm(E_far):.4e}")
# 应按 1/r 衰减
# 5. 介质界面处切向E连续性
# 在介质表面两侧取点，检查切向分量差
for test_point_pair in interface_test_points:
r_in, r_out = test_point_pair
E_in = compute_total_field(r_in, E_sca)
E_out = compute_total_field(r_out, E_sca)
n_hat = get_surface_normal(r_in)
E_tan_in = E_in - np.dot(E_in, n_hat) * n_hat
E_tan_out = E_out - np.dot(E_out, n_hat) * n_hat
jump = np.linalg.norm(E_tan_in - E_tan_out)
print(f"界面切向跳跃: {jump:.4e}") # 应趋于0
# 6. 各向异性效应：检查交叉极化
rcs_co = compute_rcs_component(theta=np.pi, phi=0, pol='xx') # co-pol
rcs_cross = compute_rcs_component(theta=np.pi, phi=0, pol='yx') # cross-pol
print(f"后向 co-pol RCS: {rcs_co:.2f} dBsm")
print(f"后向 cross-pol RCS: {rcs_cross:.2f} dBsm")
# 各向异性时 cross-pol 应该非零
# 各向同性时 cross-pol 应该极小 (< -40 dB)
二十三、完整的公式速查卡
23.1 一览表：$e^{-j\omega t}$，各向异性无磁介质，四面体一阶棱边元
┌─────────────────────────────────────────────────────────────┐
│ 完 整 公 式 速 查 卡 │
│ e^{-jωt}，各向异性无磁介质，四面体棱边元 │
├─────────────────────────────────────────────────────────────┤
│ │
│ 【波动方程】 │
│ ∇×(∇×E^sca) - k₀²ε̄̄_r E^sca = k₀²(ε̄̄_r - Ī)E^inc │
│ │
│ 【弱形式】 │
│ ∫(∇×T)·(∇×E)dV - k₀²∫T·ε̄̄_r·E dV │
│ - jk₀∮_{ABC}(n̂×T)·(n̂×E)dS │
│ = k₀²∫_{Ωd} T·(ε̄̄_r-Ī)E^inc dV │
│ │
│ 【ABC（e^{-jωt}）】 │
│ n̂×(∇×E) = jk₀ n̂×(n̂×E) 在 S_{ABC}上 │
│ │
│ 【棱边基函数】（节点 i,j，i<j） │
│ N_{ij} = λ_i∇λ_j - λ_j∇λ_i │
│ ∇×N_{ij} = 2∇λ_i×∇λ_j （常矢量） │
│ │
│ 【体积坐标梯度】 │
│ [J] = [r₁-r₄; r₂-r₄; r₃-r₄] │
│ (∇λ₁,∇λ₂,∇λ₃)ᵀ = J^{-T} │
│ ∇λ₄ = -(∇λ₁+∇λ₂+∇λ₃) │
│ V^e = |det J|/6 │
│ │
│ 【刚度矩阵】 │
│ S^e_{ij} = 4(∇λ_a×∇λ_b)·(∇λ_c×∇λ_d)·V^e │
│ 棱边 i↔(a,b)，棱边 j↔(c,d) │
│ │
│ 【各向异性质量矩阵】 │
│ T^e_{ij} = D_{bibj}I_{aiaj} - D_{biaj}I_{aibj} │
│ - D_{aibj}I_{biaj} + D_{aiaj}I_{bibj} │
│ │
│ D_{pq} = (∇λ_p)ᵀ ε̄̄_r (∇λ_q) │
│ I_{pq} = (1+δ_{pq})V^e/20 │
│ │
│ 【ABC面矩阵】 │
│ B^f_{αβ} = -jk₀[同形式展开，用 J^f_{mn}] │
│ J^f_{mn} = (1+δ_{mn})A^f/12 │
│ 注：面上节点的λ在面上非零，对面节点λ=0 │
│ │
│ 【组装】 │
│ A_{I,J} += s_i·s_j·(S^e_{ij} - k₀²T^e_{ij}) │
│ b_I += s_i·f^e_i │
│ │
│ s_k = +1 若局部方向与全局一致 │
│ s_k = -1 若方向相反 │
│ 全局棱边：(min(N_A,N_B), max(N_A,N_B)) │
│ │
│ 【入射场（e^{-jωt}，+z传播，x极化）】 │
│ E^inc = E₀ x̂ e^{+jk₀z} │
│ │
│ 【源项】 │
│ f^e_i = k₀²∫_{Ωe} N_i·(ε̄̄_r-Ī)E^inc dV │
│ 仅在介质单元内非零 │
│ │
│ 【积分公式】 │
│ 四面体：∫λ₁^a λ₂^b λ₃^c λ₄^d dV │
│ = a!b!c!d!/(a+b+c+d+3)! · 6V^e │
│ 三角形：∫λ₁^a λ₂^b λ₃^c dS │
│ = a!b!c!/(a+b+c+2)! · 2A^f │
│ │
│ 【场恢复】 │
│ E^{(e)}_k(局部) = s_k · E_{I(k)}(全局) │
│ E(r) = Σ_{k=1}^{6} E^{(e)}_k N_k(r) │
│ E^tot = E^sca + E^inc │
│ │
│ 【局部棱边定义】 │
│ e0=(0,1) e1=(0,2) e2=(0,3) │
│ e3=(1,2) e4=(1,3) e5=(2,3) │
│ │
│ 【局部面定义】（对面节点，面上节点） │
│ f0: 对面n0, 面(n1,n2,n3), 棱e3,e4,e5 │
│ f1: 对面n1, 面(n0,n2,n3), 棱e1,e2,e5 │
│ f2: 对面n2, 面(n0,n1,n3), 棱e0,e2,e4 │
│ f3: 对面n3, 面(n0,n1,n2), 棱e0,e1,e3 │
│ │
└─────────────────────────────────────────────────────────────┘
23.2 常见错误检查表
错误 症状 检查方法 忘记方向符号 $s_i s_j$ 解不收敛或完全错误 检查相邻单元共享棱边的自由度是否一致 $e^{-j\omega t}$ vs $e^{+j\omega t}$ 混用 ABC产生增长模而非衰减 检查ABC面上场是否向外衰减 $\nabla\lambda_4$ 计算错误 基函数不满足分片常数条件 验证 $\sum_i \nabla\lambda_i = 0$ Jacobian符号（正负体积） 法向量方向错误 确保 $\det J > 0$ 或统一处理 ABC面法向量朝内 能量增长而非吸收 检查法向量指向远离散射体方向 介质界面未对齐网格 虚假散射 确保介质边界与单元面重合 各向同性当各向异性处理 无误但效率低 自由空间可用简化公式 质量矩阵 $I_{pq}$ 公式混淆体/面 数值错误 体：$(1+\delta)/20\cdot V$；面：$(1+\delta)/12\cdot A$ edge_map 和 sign_map 两个数组完整描述