混合装配¶

许多偏微分方程系统在一个弱形式中耦合多个未知场：斯托克斯/纳维-斯托克斯方程中的速度和压力，孔隙弹性中的位移和压力，布辛涅斯克对流中的速度、压力和温度。对它们进行稳定离散化通常要求每个场使用不同的函数空间——经典例子是泰勒-胡德元对（二次速度、线性压力），它满足等阶元对违反的 inf-sup（LBB）条件。

MixedElementAssembler 将此类多场双线性形式装配为一个分块稀疏矩阵。你只需声明一次场，将耦合弱形式的标量被积函数写成一个 forward，后续所有操作——通过 Condenser 施加边界条件、.solve()、autograd——都在分块系统上工作，与单场情形完全一致：

from tensormesh import Condenser, Field, Mesh, MixedElementAssembler

class StokesAssembler(MixedElementAssembler):
    fields = [
        Field(trial="u", test="v", order=2, components=2),  # P2 velocity
        Field(trial="p", test="q", order=1),                # P1 pressure
    ]

    def __post_init__(self, mu=1.0):
        self.mu = mu

    def forward(self, gradu, p, gradv, q):
        return self.mu * (gradu * gradv).sum() \
            - p * gradv.diagonal().sum() \
            - q * gradu.diagonal().sum()

mesh = Mesh.gen_rectangle(chara_length=0.05, order=2).double()
assembler = StokesAssembler.from_mesh(mesh, mu=1.0)
K = assembler()          # SparseMatrix over all velocity + pressure DOFs

本页余下部分逐步介绍各个组成部分：场声明、被积函数约定、分块自由度布局、载荷向量、数据通道，以及广义阶数对——多项式阶数与网格阶数不同的场。

声明场¶

fields 是一个类属性：Field 声明的列表，每个未知场一个。每个场命名其试函数和检验函数（即你的 forward 将使用的参数名）、多项式 order，以及向量 components 的数量：

fields = [
    Field(trial="u", test="v", order=2, components=2),
    Field(trial="p", test="q", order=1),
]

三条需要记住的约定：

试函数索引列，检验函数索引行。 分块 \(K_{\beta\alpha}\) 将试场 \(\alpha\)（列）耦合到检验场 \(\beta\)（行）。对于对称弱形式你不会注意到这一点；对于非对称耦合（对流、散度约束），该声明显式地确定了方向。
名称是全局的。 所有试函数和检验函数名称必须是互异的有效标识符，不能为 "x"（坐标保留字），也不能以 "grad" 开头（梯度派生为 "grad" + name）。
阶数是按场的，与网格阶数无关——参见广义阶次配对。

标量被积函数¶

forward 返回耦合双线性形式在单个求积点处的标量被积函数——与 ElementAssembler 的按名分派相同（参见弱形式约定），但场参数现在是张量值的：

参数	每点形状	含义
标量场值（`p`、`q`）	0维 `[]`	求积点处的场值
标量场梯度（`gradp`）	`[D]`	物理空间梯度
向量场值（`u`、`v`）	`[c]`	全部 \(c\) 个分量一起
向量场梯度（`gradu`）	`[c, D]`	雅可比矩阵：第 \(i\) 行是 \(\nabla u_i\)

因此 gradu.diagonal().sum() 是 \(\nabla\cdot u\)，(gradu * gradv).sum() 是 \(\nabla u : \nabla v\)，而 (gradu @ w).dot(v) 是对流项 \((w\cdot\nabla)u\cdot v\)。仅计算你命名的参数：若某个分块的传递不涉及任何场，则跳过该分块，因此缺失的耦合（例如斯托克斯方程的压力-压力零分块）除了稀疏模式中的占位外不产生任何开销。

重要

被积函数必须是双线性的：每一项必须恰好包含一个试函数因子和一个检验函数因子。内部通过对其他场置零来提取每个 \((\alpha, \beta)\) 分块——这对双线性形式是恒等操作，否则则是静默错误。装配器通过廉价的零求值检查来防范这一点，若被积函数含常数项或非双线性部分则抛出异常。非线性弱形式按有限元常规方式处理：线性化（皮卡/牛顿）并将前一次迭代作为数据传入，例如下面滞后的对流速度 w。

例如，皮卡线性化的纳维-斯托克斯步在斯托克斯被积函数上增加一个数据参数（完整脚本参见顶盖驱动方腔）：

def forward(self, gradu, p, v, gradv, q, w):
    convection = self.rho * (gradu @ w).dot(v)
    diffusion = self.mu * (gradu * gradv).sum()
    return convection + diffusion \
        - p * gradv.diagonal().sum() \
        - q * gradu.diagonal().sum()

构建与调用¶

from_mesh() 构建装配器；额外的位置/关键字参数传给你的 __post_init__（物理参数），与单场装配器完全一致：

assembler = StokesAssembler.from_mesh(mesh, mu=1.0)
K = assembler()                       # -> SparseMatrix  [N, N]

默认求积阶数为 2 * max(field orders)——对仿射单元上两个最高阶基函数的乘积精确。传入 quadrature_order= 覆盖（例如用于数据加权项或对流项）；内置表最高到7阶。

结果是一个大小为 \(N = \sum_f n_f\, c_f\) 的方阵 SparseMatrix。重复调用（皮卡迭代、时间步）复用相同的缓存索引张量，因此一次构建的 Condenser 在多次调用间保持其凝聚布局不变。

备注

混合装配目前为单设备：distributed() 装饰器（分布式有限元方法）仅覆盖单场装配器。跨进程分布分块布局需要尚不存在的全局自由度编号层——这是一项计划中的 ROADMAP 条目。

分块自由度布局¶

全局自由度向量按声明顺序堆叠各场，每场内按节点优先：

\[\mathrm{dof} \;=\; \mathrm{offset}_f + n_{\mathrm{local}}\cdot c_f + \mathrm{comp}.\]

你无需手动索引——assembler.layout（一个 BlockLayout）提供了相关操作：

layout = assembler.layout

layout.n_dofs                      # total N
layout.names                       # ["u", "p"]
layout.n_nodes("u")                # nodes of the P2 velocity space
layout.points("u")                 # their coordinates  [n_u, D]

# masks and indices (all in global block-DOF numbering)
layout.dof_mask("u")                          # every u DOF
layout.dof_mask("u", node_mask, component=0)  # u_x on selected nodes
layout.dof_mask("p", where=lambda x: x[:, 0] < 1e-6)   # coordinate predicate
layout.boundary_mask("u")                     # topological boundary DOFs
layout.dof_index("p", int(layout.node_ids("p")[0]))    # a single DOF (pressure pin)

# moving between the block vector and per-field tensors
fields = layout.split(sol)         # {"u": [n_u, 2], "p": [n_p]}
sol = layout.cat(u=vel, p=0.0)     # assemble a block vector
p_on_mesh = layout.prolong("p", fields["p"])   # P1 -> all mesh points
f_on_p = layout.restrict("p", point_values)    # mesh points -> P1 nodes

边界条件则通过不变的 Condenser 施加——混合系统只是一个更大的带布尔掩码的自由度向量：

bc_mask = layout.dof_mask("u", mesh.boundary_mask)     # no-slip walls
bc_mask[layout.dof_index("p", int(layout.node_ids("p")[0]))] = True

bc_val = torch.zeros(layout.n_dofs, dtype=torch.float64)
bc_val[layout.dof_mask("u", is_lid, component=0)] = 1.0

condenser = Condenser(bc_mask, bc_val[bc_mask])
K_, f_ = condenser(K, f)
sol = condenser.recover(K_.solve(f_))

boundary_mask 通过拓扑方式检测边界（仅被一个单元引用的面），因此它也适用于没有 is_boundary 点数据的网格——例如 gmsh 导入的网格——以及自由度不是网格节点的广义阶次场。

载荷向量¶

assemble_vector() 是混合空间的线性形式 \(l(v)\)：相同的分发机制，仅测试参数，输出相同块布局的稠密 [N] 向量。为问题的"自然"右端项重写 forward_vector，或传入一次性 func=：

class TransientNS(MixedElementAssembler):
    fields = [...]

    def forward(self, u, gradu, p, v, gradv, q, w):
        ...  # backward-Euler matrix side

    def forward_vector(self, v, uprev):
        return self.rho / self.dt * uprev.dot(v)

f = assembler.assemble_vector(field_data={"uprev": ("u", u_prev)})

# one-off linear forms: e.g. L2-project the vorticity of a P2
# velocity onto the P1 pressure space
omega_rhs = assembler.assemble_vector(
    func=lambda q, gradw: (gradw[1, 0] - gradw[0, 1]) * q,
    field_data={"w": ("u", velocity)},
)

由于载荷在每个场自身的空间中积分，广义阶次场上的源项按构造就是协调的——无需在之间手动插值。

数据通道¶

所有单场通道在 forward / forward_vector 中均可用——x、point_data（以网格阶次基函数插值，grad{key} 可用）、element_data、scalar_data——外加一个新通道：

field_data={"name": (field, values)}

数据位于场自身的自由度上，以该场的基函数插值（grad{name} 也可用）。这是用于 Picard / 时间步进在广义阶次场上的通道：拓扑 P2 速度的前一次迭代没有网格点表示，因此必须承载于自身的自由度向量：

w = layout.split(sol)["u"]                    # [n_u, c] on the P2 DOFs
K = assembler(field_data={"w": ("u", w)})

当网格阶次等于场阶次时，两个通道一致——此时 point_data 是更简单的写法。

广义阶次配对¶

场阶次不绑定于网格阶次。除阶次 1 和网格阶次（其自由度位于网格节点）外，任何其他阶次都会构建一个拓扑拉格朗日自由度映射：每个顶点一个自由度，每条唯一（定向）边 \(k-1\) 个，加上单元内部。实际意味着：

线性 gmsh 网格上的 Taylor-Hood。 在导入的 P1 网格上，Field(order=2) 拓扑地创建二次空间——无需在阶次 2 重新划分网格。圆柱绕流和障碍物绕流示例正是以这种方式运行。
更高阶配对。 P3-P2 及更高阶同理；边自由度存储在定向边上，因此相邻单元一致，空间在任意阶次保持 \(C^0\) 协调。
几何始终保持网格阶次的等参数。 P2 场在弯曲（阶次 2）网格上积分时即使与 P1 压力配对也使用弯曲几何；次参数和超参数场都是正确的。

当前限制，以显式错误检查：三维中拓扑映射覆盖顶点 + 边自由度（四面体 P2 可用；需要面自由度的阶次会抛出 NotImplementedError），求积表止于 7 次——实际约 P3-P2。

对于拓扑自由度场，layout.node_ids（网格节点概念）不可用；请使用 layout.points(name)、layout.boundary_mask(name) 和 layout.dof_mask(name, where=...)，它们适用于每种场类型，并用 layout.prolong / layout.restrict（等参数插值）传输数据。

实例教程¶

流体示例库完全基于混合装配构建，按复杂度递增：

Taylor-Hood Stokes — 带收敛性研究的 manufactured Stokes 解；最佳入门点。
顶盖驱动方腔 — 定常 Navier-Stokes（Picard），二维和三维。
圆柱绕流（涡脱落） — P1 gmsh 网格上使用 forward_vector 和 field_data 的瞬态 Navier-Stokes。
瑞利-贝纳德对流 — 三场 Boussinesq 系统（速度、压力、温度），含非对角浮力块。