应用主从节点法施加多点约束
有限元离散后，要求节点处内力与外力平衡，\(\boldsymbol{\mathsf{f}}^{\mathrm{int}}=\boldsymbol{\mathsf{f}}^{\mathrm{ext}}\)。其中，内力可由刚度矩阵 \(\boldsymbol{\mathsf{K}}\) 和节点自由度列阵 \(\boldsymbol{\mathsf{d}}\) 给出，最终得到离散后的平衡方程： \[ \begin{equation} \boldsymbol{\mathsf{K}} \boldsymbol{\mathsf{d}} = \boldsymbol{\mathsf{f}}^{\mathrm{ext}}. \end{equation} \] 除平衡方程之外，节点之间一般还需满足若干运动学约束： \[ F_{k}( \boldsymbol{d}_{a}, \ldots ) = \boldsymbol{0}, \quad k=1,\ldots \] 当上述运动学约束只与单个节点 \(\boldsymbol{d}_{a}\) 的自由度相关时，此时等价于在该节点处施加位移边界条件。当运动学约束与多个节点相关时，此时称该约束为多点约束。由于这些运动学约束的存在，刚度矩阵 \(\boldsymbol{\mathsf{K}}\) 一般不是满秩矩阵，需要在施加约束之后才能唯一定解节点位移。施加位移边界条件 (单点约束) 可直接消除刚度矩阵对应自由度编号的行和列，而施加多点约束一般有如下三种方式：(1) 主从节点法 (2) 罚函数法 和 (3) 拉格朗日乘子法。
主从节点法是将多点约束直接代入到线性方程组中，消除从节点自由度。虽然该方法思路简单，但由于
- 确定主从节点对的方案并不唯一，并且很难找到一个通用的算法，在任意给定的约束方程组中确定正确约束数量的主从节点对;
- 施加约束需要对刚度矩阵的元素进行调整，移动或删除操作，这可能会在刚度矩阵按照稀疏矩阵存储时遇到困难;
- 该方法只能施加线性约束条件,
所以主从节点法在通用性方面远不如 罚函数法 或 拉格朗日乘子法。不过在针对特定问题时 (例如，求解在周期性边界条件约束下单胞内的线弹性问题)，多点约束的形式较单一，所以仍可以找到一套系统的算法，确定问题中的主从节点对和高效调整刚度矩阵。
主从节点法
主从节点法本质上是代入方法，对于原始未知量 \(\boldsymbol{\mathsf{d}}\)，维度为
NDOF x 1，在施加 NcDOF 约束之后，可通过数量为
NiDOF=NDOF-NcDOF 的独立未知量 \(\hat{\boldsymbol{\mathsf{d}}}\) 表示为
\[
\begin{equation}
\underbrace{\boldsymbol{\mathsf{d}}}_{\scriptsize \mathrm{NDOF}\times1}
= \underbrace{\boldsymbol{\mathsf{T}}}_{\scriptsize
\mathrm{NDOF}\times\mathrm{NiDOF}}\
\underbrace{\hat{\boldsymbol{\mathsf{d}}}}_{\scriptsize
\mathrm{NiDOF}\times1}.
\end{equation}
\] 将该方程代入平衡方程之后，再左乘 \(\boldsymbol{\mathsf{T}}^{\top}\)，就得到满秩的线性代数方程组：
\[
\begin{equation}
\underbrace{\hat{\boldsymbol{\mathsf{K}}}}_{\scriptsize
\mathrm{NiDOF}\times\mathrm{NiDOF}}\
\underbrace{\hat{\boldsymbol{\mathsf{d}}}}_{\scriptsize
\mathrm{NiDOF}\times1} =
\underbrace{\hat{\boldsymbol{\mathsf{f}}}}_{\scriptsize
\mathrm{NiDOF}\times1}, \quad \text{式中}\quad
\hat{\boldsymbol{\mathsf{K}}} = \boldsymbol{\mathsf{T}}^{\top}
\boldsymbol{\mathsf{K}} \boldsymbol{\mathsf{T}}, \quad
\hat{\boldsymbol{\mathsf{f}}} = \boldsymbol{\mathsf{T}}^{\top}
\boldsymbol{\mathsf{f}}.
\end{equation}
\]
约束自由度又可以分成两类：施加位移边界条件的自由度 和
多点约束中的从节点自由度，其数量分别为 NdcDOF 和
NscDOF。因此，最终方程中独立未知数的个数为 \[
\mathrm{NiDOF} = \mathrm{NDOF} - \mathrm{NdcDOF} - \mathrm{NscDOF}
\]
使用主从节点法的主要难点在于：
- 确保约束矩阵 \(\boldsymbol{\mathsf{T}}\) 是满秩矩阵。
- 在代码实现中避免构造 \(\boldsymbol{\mathsf{T}}\)，获取修正后的刚度矩阵 \(\hat{\boldsymbol{\mathsf{K}}}\)。
以下将针对线下阶段影响张量的求解，分别给出对上述两个问题的具体实现。
应用：线下阶段求解影响张量
在对方形单胞相对面 \(S^{-}\)，\(S^{+}\) 施加周期性边界条件时，方形单胞的相对面节点位移满足如下约束： \[ \boldsymbol{d}_{a}^{+}-\boldsymbol{d}_{a}^{-} = \boldsymbol{d}_{b}^{+} - \boldsymbol{d}_{b}^{-}, \quad \forall ~\text{node}~ a,b \in S \]
单胞问题主从节点对的确定
在求解影响张量时，施加在单胞上的宏观应变等于零，这意味着对于方形单胞，8
个角点的位移全部等于零： \[
\boldsymbol{d}_{a}=0, \quad \text{node}\ a \text{ at corner}
\] 正是因为宏观应变等于零，所以对任意相对面上的节点，位移均相等：
\[
\begin{equation}
\boldsymbol{d}_{a}^{+} = \boldsymbol{d}_{a}^{-}, \quad \forall
~\text{node}~ a,b \in S
\end{equation}
\] 因此，确定约束矩阵 \(\boldsymbol{\mathsf{T}}\)，就转变为寻找所有在
(1) 边界面 和 (2) 棱线 处的节点对 \(\{ a^{+},
a^{-} \}\)，\(\{a^{++}, a^{+-}, a^{-+},
a^{--}\}\)，并在节点对中确定一个节点作为主节点，其它所有节点作为从节点，将从节点的三个自由度全部设置为属性
slave，将在之后调整线性方程组时消去该自由度。
线性方程组的调整
首先，对于属于位移边界条件的自由度，在组装全局刚度矩阵时直接跳过，由此得到数量为
NDOF-NdcDOF 的线性方程组： \[
\boldsymbol{\mathsf{K}}_{uu} \boldsymbol{\mathsf{d}}_{u} =
\boldsymbol{\mathsf{f}}_{u}
\] 将未知量 \(\boldsymbol{\mathsf{d}}_{u}\) 分解为 \((\hat{\boldsymbol{\mathsf{d}}};
\boldsymbol{\mathsf{d}}_{s})\)，其中 \(\hat{\boldsymbol{\mathsf{d}}}\)
表示独立的自由度，维度为 NiDOF x 1，\(\boldsymbol{\mathsf{d}}_{s}\)
表示从属自由度，维度为 NscDOF x 1。对于式 (4)
中的约束关系，可以定义维度为 NscDOF x NiDOF 矩阵 \(\boldsymbol{\mathsf{A}}\)，在第
sidx 行处只在列索引为 s2mMap(sidx)
处等于1，其它位置均等于零，s2mMap
将从属自由度索引映射到对应的主自由度索引。由此，式 (4)
约束关系，以及约束矩阵 \(\boldsymbol{\mathsf{T}}\) 可以表示为 \[
\boldsymbol{\mathsf{d}}_{s} = \boldsymbol{\mathsf{A}}
\hat{\boldsymbol{\mathsf{d}}}, \quad
\boldsymbol{\mathsf{T}} = \begin{bmatrix}\boldsymbol{\mathsf{I}} \\
\boldsymbol{\mathsf{A}} \end{bmatrix}
\] 所以，添加约束关系后的刚度矩阵和外力列阵等于 \[
\hat{\boldsymbol{\mathsf{K}}} = \begin{bmatrix}\boldsymbol{\mathsf{I}}
& \boldsymbol{\mathsf{A}} ^{\top}\end{bmatrix}
\boldsymbol{\mathsf{K}} \begin{bmatrix}\boldsymbol{\mathsf{I}} \\
\boldsymbol{\mathsf{A}} \end{bmatrix}, \quad
\hat{\boldsymbol{\mathsf{f}}} = \begin{bmatrix}\boldsymbol{\mathsf{I}}
& \boldsymbol{\mathsf{A}} ^{\top}\end{bmatrix}
\boldsymbol{\mathsf{f}}.
\] 此时，\(\boldsymbol{\mathsf{A}}^{\top}\)
的作用相当于将索引为 sidx 的从属自由度组装在
s2mMap(sidx) 位置。