静磁学（麦克斯韦方程）¶
一个三维静磁学问题：通过在单位立方体上求解磁 矢量势 \(\mathbf{A}\)，恢复载流导线周围的磁场。脚本 examples/maxwell/magnetostatic.py 采用连续节点拉格朗日单元，并使用 Badia 与 Codina（2012）[BC2012] 的 稳定化混合格式，使旋度-旋度问题在标准的（非旋度协调的）节点空间上适定。
备注
TensorMesh 自带连续的 节点拉格朗日 单元。\(\nabla\times\nabla\times\) 的自然空间是 旋度协调 的 Nédélec 边单元，目前尚未原生支持——因此本示例采取了与 顶盖驱动方腔 相同的务实路线：在等阶节点空间上求解，并添加 一致稳定化（这里是库仑规范的拉格朗日乘子，加上 grad-div 项和乘子项）以抑制伪模。更专门的有限元空间将逐步加入。
问题¶
静磁学由安培定律 \(\nabla\times\mathbf{H}=\mathbf{J}\)，连同 \(\nabla\cdot\mathbf{B}=0\) 和 \(\mathbf{B}=\mu\mathbf{H}\) 共同描述。写成 \(\mathbf{B}=\nabla\times\mathbf{A}\) 可恒等地满足散度约束，并且——将磁阻率 \(\nu=1/\mu\) 吸收进一个无量纲缩放后——将安培定律约化为关于 \(\mathbf{A}\) 的 旋度-旋度 方程：
其中标量 \(p\) 是一个拉格朗日乘子，用于施加 库仑规范 \(\nabla\cdot\mathbf{A}=0\)。该规范并非可有可无：旋度-旋度算子会湮灭任意梯度场，因此其核空间极其庞大，没有规范时 \(\mathbf{A}\) 仅在相差一个 \(\nabla\phi\) 的意义下唯一。
边界条件为齐次切向迹
即在每个与坐标轴对齐的面上，\(\mathbf{A}\) 的两个 切向 分量被固定为零（在 \(x=0,1\) 上即 \(A_y\) 和 \(A_z\)），这些约束沿边和角处叠加组合。
电流源。 一条穿过立方体中心、沿 \(z\) 方向的光滑通道——直导线的一个柔化模型，
其中 \(J_0 = 100\)，中心 \((x_c, y_c) = (0.5, 0.5)\)，宽度 \(\sigma = 0.08\)。它沿 \(z\) 方向不变，因此 \(\nabla\cdot\mathbf{J} = 0\)；在自由空间中，这样的电流会产生一个纯 方位角 方向、环绕导线的场，这正是我们要寻找的定性特征。
稳定化弱形式¶
在连续拉格朗日空间中求 \((\mathbf{A}_h, p_h) \in V_h \times Q_h\)，使得对所有 \((\mathbf{v}_h, q_h)\) 均有
其中各双线性型为
正是这两个稳定化项使得节点离散稳健：grad-div 惩罚项 \(s_u\) 控制 \(\mathbf{A}\) 的散度，而 \(s_p\) 对乘子进行正则化——这就是 [BC2012] 的方程 (3.6)。该演示对每个单元使用常值长度尺度 \(h^2 = \texttt{chara\_length}^2\)。装配之后，这四个形式构成一个 2×2 分块鞍点系统
TensorMesh 设置¶
每个双线性型都是一个小型的 ElementAssembler，其 forward 返回逐节点的分块。例如，旋度-旋度形式构造两个反对称的叉乘矩阵 \([\nabla N]_\times\) 并对它们做缩并——\([\nabla N]_\times\) 的第 \(a\) 列恰好是 \(\nabla\times(N\,\mathbf{e}_a)\)，因此该缩并就是 \([3,3]\) 分块 \((\nabla\times\mathbf{u})\cdot(\nabla\times\mathbf{v})\)：
class CurlCurlAssembler(ElementAssembler):
"""(curl u, curl v) for 3D vector nodal fields."""
def forward(self, gradu, gradv):
zero = torch.zeros_like(gradu[0])
curl_u = torch.stack((torch.stack((zero, -gradu[2], gradu[1])),
torch.stack((gradu[2], zero, -gradu[0])),
torch.stack((-gradu[1], gradu[0], zero))))
curl_v = ... # the same, built from gradv
return curl_u.T @ curl_v # [3, 3] block
A = CurlCurlAssembler.from_mesh(mesh)()
Su = DivergenceStabilizationAssembler.from_mesh(mesh, h2=h2)() # h^2 (div u, div v)
Sp = PressureStabilizationAssembler.from_mesh(mesh)() # (grad p, grad q)
B = PressureCouplingAssembler.from_mesh(mesh)() # -(grad p, v)
K = SparseMatrix.combine([[A + Su, B],
[-1.0 * B.T, Sp]])
几处值得指出的细节：
分块装配。
SparseMatrix.combine将四个稀疏分块拼贴成鞍点矩阵，使得矢量 \(\mathbf{A}\) 自由度和标量 \(p\) 自由度共享同一个线性系统。作为自由度掩码的切向边界条件。
tangential_vector_potential_mask为 \(\mathbf{n}\times\mathbf{A}=\mathbf{0}\) 构建逐分量的狄利克雷掩码；将其与用于 \(p\) 的mesh.boundary_mask拼接后，喂入单个Condenser，在求解前凝聚掉每一个受约束的自由度。恢复场。 \(\mathbf{B}=\nabla\times\mathbf{A}\) 是一个导出的、逐单元的量。脚本通过 \(L^2\) 投影将其投回到光滑的节点场——装配一个
MassElementAssembler，用一个NodeAssembler构建旋度右端项，再逐分量求解 \(M\,\mathbf{B} = \tilde{\mathbf{B}}\)。这与 圆柱绕流 恢复涡量所用的方法相同。
图 50 magnetostatic.py 的输出（在 ParaView 中渲染）：\(\mathbf{B}=\nabla\times\mathbf{A}\) 的流线，大致沿 \(z\) 轴向下观察。中心处暖色的图元标示出沿 \(z\) 方向的电流通道；场线沿方位角方向环绕它——直导线的右手定则图样——而有限的立方体边界则将外圈轻微地方形化。¶
运行它¶
cd examples/maxwell
python magnetostatic.py # writes magnetostatic_3d.vtu
求解会写出 magnetostatic_3d.vtu，其中包含两个节点点场：矢量势 A 和磁场 curl_A \(= \nabla\times\mathbf{A}\)。在 ParaView 中打开它，并对 curl_A 应用 Stream Tracer 或 Glyph 即可复现上图。如果你想做实验，网格大小和电流设置就在 main() 的开头。