边界条件¶

当你从形式得到了刚度矩阵 K 和载荷向量 b 之后，仍需施加问题的边界条件。在如今的 TensorMesh 中，这是通过对狄利克雷自由度做 静态凝聚 实现的——只需对 Condenser 调用一次，便能接收一个完整系统，并返回一个仅在内部自由度上的、可直接求解的缩减系统。

诺伊曼条件不需要单独的算子：它们作为边界积分自然地出现在弱形式中，由 FacetAssembler 装配。

本章其余部分将逐一讲解这两者，外加向量值问题（线弹性、Stokes 等）所需的逐分量掩码、随时间变化的边界数据，以及实践中最常遇到的陷阱。

为何采用静态凝聚¶

将自由度划分为“内部”（自由）和“外部”（被狄利克雷固定）后，可将线性系统

\[\begin{split}\begin{bmatrix} K_{ii} & K_{io} \\ K_{oi} & K_{oo} \end{bmatrix} \begin{bmatrix} u_i \\ u_o \end{bmatrix} = \begin{bmatrix} b_i \\ b_o \end{bmatrix}\end{split}\]

转化为仅在内部块上的一个更小求解：

\[K_{ii} \, u_i \;=\; b_i - K_{io} \, u_o, \qquad u_o \text{ prescribed.}\]

凝聚避免了拉格朗日乘子，能让对称正定（SPD）算子保持 SPD 性质，并且端到端可微——边界值会经由右端项的修正流动，并以通常的方式获得梯度。

Condenser API¶

从一个覆盖所有自由度的布尔掩码构造 Condenser：

from tensormesh import Condenser

condenser = Condenser(mesh.boundary_mask)              # zero values (default)
condenser = Condenser(mesh.boundary_mask, values)      # prescribed values

dirichlet_mask: [n_dof] —— 布尔张量；在自由度被固定处取 True。对于标量问题，它就是 mesh.boundary_mask（长度为 n_points）；对于向量值问题，参见下文的向量值问题（弹性、Stokes）。
dirichlet_value —— 可选的一维张量，给出指定的边界值。它可以是 [n_dof]（构造函数会将其切片至边界项）或已切片好的 [n_outer_dof]。默认为全零。

三个方法即可覆盖常见的工作流程：

方法	作用
`condenser(K, b)`	一次性凝聚两者。返回 `(K_inner, b_inner)`。请先运行它——它会为后续调用缓存稀疏布局。
`condenser.condense_rhs(b_new)`	仅重新凝聚右端项，复用缓存的布局。当 `K` 固定而 `b` 变化时（瞬态问题、批量右端项）很有用。
`condenser.recover(u_inner)`	将内部解与指定的边界值重新拼接成一个完整的 `[n_dof]` 向量。
`condenser.restrict(f)` / `condenser.prolong(f_inner)`	不涉及边界值的线性投影（在边界上以零进行丢弃/散布）。当所提升的量无论狄利克雷值如何都应在边界上为零时使用——例如 `ImplicitLinearRungeKutta` 中各阶段的斜率。积分器类的相关模式参见与边界条件组合。

对于含时边界条件，condenser.update_dirichlet(new_values) 会换入一个新的边界值向量（与构造函数一样，既接受完整的 [n_dof] 形式，也接受已切片好的 [n_outer_dof] 形式）。缓存的稀疏布局与所指定的值无关，因此该调用开销很小。

首次调用 condenser(K, b) 会计算并缓存 K 稀疏模式的内部/外部索引划分；后续对具有相同布局的矩阵的调用会复用这一成果。

备注

Condenser 是一个 torch.nn.Module，它所拥有的每个张量（布尔掩码、指定值，以及惰性计算的索引缓冲区）都注册为缓冲区。因此

condenser = Condenser(mesh.boundary_mask).to(mesh.device)

可将该算子与 K、b 和网格一起迁移。如果你忘了这么做，__call__ 和 condense_rhs() 仍会在每次调用时自动提升到输入矩阵的设备/数据类型——这很方便，但在紧凑的瞬态循环中，逐次调用的 .to(...) 开销会在性能分析中显现出来。

齐次狄利克雷¶

最基本、最常见的情形——边界上取零值——正是快速开始所做的：

condenser = Condenser(mesh.boundary_mask)
K_, b_ = condenser(K, b)
u_inner = K_.solve(b_)
u = condenser.recover(u_inner)

连同求解在内只需三行。恢复出的 u 在边界自由度处为零，在其余各处为有限元解。

非齐次狄利克雷¶

通过将逐自由度的值传给构造函数来指定它们。可接受两种形状——选用更方便的一种即可：

# Option A: full-length [n_dof] vector (Condenser slices internally)
values = torch.zeros(mesh.n_points, dtype=mesh.dtype)
x = mesh.points[:, 0]
values[(x == 1.0) & mesh.boundary_mask] = 1.0       # u = 1 on right edge
condenser = Condenser(mesh.boundary_mask, values)

# Option B: pre-sliced [n_outer_dof] vector
sliced = values[mesh.boundary_mask]
condenser = Condenser(mesh.boundary_mask, sliced)

凝聚后的右端项会自动吸收项 \(-K_{io}\, u_o\)，因此流程的其余部分保持不变。

小技巧

Mesh.gen_rectangle / gen_cube 将边界节点恰好放置在整数坐标处，因此 x == 1.0 是有效的。对于从 Gmsh / VTK 文件加载的网格，建议改用

mask = torch.isclose(x, torch.tensor(1.0, dtype=mesh.dtype)) & mesh.boundary_mask

这样网格生成器带来的 1e-16 量级扰动就不会在不报错的情况下丢掉一排节点。

向量值问题（弹性、Stokes）¶

对于弹性问题，每个节点带有 dim 个位移分量，而由 LinearElasticityElementAssembler 产生的装配器会生成一个尺寸为 n_dof = n_points * dim 的全局系统。两个额外步骤可将 mesh.boundary_mask（形状 [n_points]）转化为自由度级别的狄利克雷掩码：

dim = mesh.points.shape[1]

# Per-(node, component) mask, then flatten to per-DOF
bc_mask = torch.zeros(mesh.n_points, dim, dtype=torch.bool)

# Clamp the bottom edge in all directions, leave top free
y = mesh.points[:, 1]
bc_mask[(y == 0.0) & mesh.boundary_mask, :] = True

# Pin x-displacement on the right edge (roller BC), y free
x = mesh.points[:, 0]
bc_mask[(x == 1.0) & mesh.boundary_mask, 0] = True

condenser = Condenser(bc_mask.flatten())          # [n_points * dim]

同样的展平技巧也适用于指定值：先构建一个 [n_points, dim] 的目标位移张量，再在交给构造函数之前 .flatten()。约定是自由度按节点交错排列（[u_x_0, u_y_0, u_x_1, u_y_1, ...]），与 TensorMesh 弹性装配器的布局相匹配。

随时间变化的边界值¶

只需构建一次凝聚器，随后便可在每个时间步推入新值而无需重建布局：

condenser = Condenser(mesh.boundary_mask, initial_values)
K_, _ = condenser(K, torch.zeros(mesh.n_points, dtype=mesh.dtype))   # caches the split

for t in time_steps:
    condenser.update_dirichlet(values_at(t))            # cheap setter
    b_new   = build_rhs_for(t)
    b_inner = condenser.condense_rhs(b_new)             # reuses cached K_io
    u       = condenser.recover(K_.solve(b_inner))

condense_rhs() 要求 __call__ 已在相同的矩阵布局上至少被调用过一次——正是这次调用填充了它所使用的缓存 \(K_{io}\) 块。请在循环之前调用它，哪怕你会丢弃第一次求解的结果。

将狄利克雷与诺伊曼结合¶

在 \(\partial\Omega\) 一部分上的非零牵引力或通量，对应弱形式右端的一个边界积分：

\[l(v) \;=\; \int_{\Omega} f \cdot v \, \mathrm{d}\Omega \;+\; \int_{\Gamma_N} g \cdot v \, \mathrm{d}S .\]

这里不涉及 Condenser——自然边界条件完全由弱形式处理。用一个 FacetAssembler 子类装配面积分，并在凝聚之前将结果加到你的载荷向量上：

import torch
from tensormesh import Mesh, Condenser, FacetAssembler
from tensormesh.assemble import LaplaceElementAssembler

mesh = Mesh.gen_rectangle(chara_length=0.1)
x    = mesh.points[:, 0]

# ----------------------- volume bilinear form ----------------------------
K = LaplaceElementAssembler.from_mesh(mesh)()

# ----------------------- Neumann flux on the right edge ------------------
right_edge = (x == 1.0) & mesh.boundary_mask         # [n_points] bool

class FluxAssembler(FacetAssembler):
    def forward(self, v):
        return v                                     # ∫_{Γ_N} g·v  with g=1

b_neumann = FluxAssembler.from_mesh(mesh, boundary_mask=right_edge)()

# Optional volume load f = 1
from tensormesh.assemble import const_node_assembler
b_volume = const_node_assembler(1.0).from_mesh(mesh)()

b = b_volume + b_neumann

# ----------------------- Dirichlet on the *complement* ------------------
dirichlet = mesh.boundary_mask & ~right_edge         # bottom/top/left edges
condenser = Condenser(dirichlet)

K_, b_ = condenser(K, b)
u      = condenser.recover(K_.solve(b_))

有两点值得强调：

FacetAssembler.from_mesh(mesh, boundary_mask=right_edge) 接受与 Condenser 所用相同类型的点级布尔掩码，因此 同一个 逻辑上的“右边缘”张量既驱动诺伊曼积分，也决定哪些节点保持自由。
传给 Condenser 的狄利克雷掩码是排除了诺伊曼区域之后的边界——否则指定值会覆盖通量积分的贡献。

混合区域与逐分量边界条件¶

对于相互独立的狄利克雷与诺伊曼区块（或弹性中的逐分量混合），可从坐标构建掩码，并用布尔代数将它们组合起来：

x, y = mesh.points[:, 0], mesh.points[:, 1]
left  = (x == 0.0) & mesh.boundary_mask
right = (x == 1.0) & mesh.boundary_mask

# Scalar problem: fix left to zero, ramp on right
dirichlet_mask = left | right
values = torch.zeros(mesh.n_points, dtype=mesh.dtype)
values[right] = 1.0
condenser = Condenser(dirichlet_mask, values)

对于从 Gmsh 加载、带有命名物理组的网格，这些组会以布尔掩码的形式落入 mesh.point_data；你可以按键读取它们（例如 mesh.point_data["clamped_face"]），而无需手动从坐标构建掩码。

易踩的坑¶

稀疏布局在首次调用时被缓存。 Condenser 会记住你首次通过 __call__ 交给它的矩阵的行/列划分；传入具有不同稀疏模式的矩阵会触发断言错误。如果你的几何或单元类型发生了变化，请新建一个 Condenser。
向量值掩码必须展平。 一个 [n_points, dim] 的布尔张量并不是有效的 dirichlet_mask——传入前请调用 .flatten()。同理，弹性问题的指定值应当是一个展平后的 [n_points * dim]（或 [n_outer_dof]）向量。
浮点扰动网格上的坐标掩码。 精确比较（x == 1.0）在 TensorMesh 的内置生成器上有效，但在导入的网格上可能失效；拿不准时请使用 torch.isclose。
重叠区域上狄利克雷优先。 如果某个节点既在狄利克雷掩码中，又落在某个诺伊曼面积分的支撑范围内，那么该行上的诺伊曼贡献会被凝聚丢弃（该行被移除）。当你希望两者都生效时，请将诺伊曼区块从狄利克雷掩码中排除。
求解与恢复必须由同一个掩码驱动。 不要在凝聚与恢复之间修改 mesh.boundary_mask——Condenser 引用的是构造时传给它的那同一个缓冲区。

目前尚未内置的功能¶

tensormesh.operator 模块恰好只提供一个算子——通过静态凝聚施加狄利克雷条件的 Condenser。还有几种相邻的技术并非一等算子：

用于狄利克雷的 拉格朗日乘子 —— 凝聚的一种替代方案，它保留原始的自由度布局，但会产生一个鞍点系统。如果你自行构建增广矩阵，手动实现是可行的，但没有现成的辅助类。
周期性边界条件 —— 通常通过识别匹配的自由度、并在装配后的矩阵中合并对应的行/列来实现。没有现成的辅助工具。
作为顶层算子的 罗宾（混合）条件 —— 不过其双线性形式贡献 \(\int_{\Gamma} \alpha\, u\, v\, \mathrm{d}S\) 只需一行代码的 FacetAssembler 子类即可；在凝聚之前将它加到 K 上就行。

罚函数接触、表面张力、压力载荷以及类似的基于能量的边界项确实是一等公民——参见 ContactAssembler（一个带有 element_energy 钩子以表示表面能密度的 FacetAssembler 子类）。示例库针对需要它们的情形提供了完整的实操方案。

接下来¶

稀疏求解器 —— 用合适的后端求解凝聚后的系统。
时间积分 —— 在瞬态问题中，同一个凝聚器会被逐步复用。
可微性 —— 对边界值本身求梯度（反问题/控制问题）。
快速开始 —— 包含齐次边界条件在内的完整流程。