这一篇讲一讲怎么将Strain作为Constraint放到XPBD的框架里。
XPBD
在前面的博客当中我已经讲过XPBD的由来和概念了。现在来回顾一下推导的步骤：
- 由隐式欧拉法得到一个运动方程
- 得到运动方程的等价形式：最小值优化问题
- 将约束和能量等价，引入拉格朗日乘子
- 得到最小化问题的解（未简化）
- 简化处理得到问题的解
运动的连续时间积分
$$\mathbf{x}^{t+1} = \mathbf{x}^{t} + \int_{t}^{t+1} \mathbf{v} dt, \\
\mathbf{v}^{t+1} = \mathbf{v}^{t} + \int_{t}^{t+1} \mathbf{a} \ dt = \mathbf{v}^{t} + \mathbf{M}^{-1} \int_{t}^{t+1} (\mathbf{f_{int}(\mathbf{x}(t)) + f_{ext}}) dt$$
运动离散时间积分隐式欧拉法
$$
\begin{align}
\mathbf{x}^{t+1} &= \mathbf{x}^{t} + h\mathbf{v}^{t+1}, \\
\mathbf{v}^{t+1} &= \mathbf{v}^{t} + h\mathbf{M}^{-1}(\mathbf{f}_{int}(\mathbf{x}_{n+1}) + \mathbf{f}_{ext})
\end{align}
$$
最小化问题
将隐式欧拉法的两个公式进行处理我们可以得到：
$$\mathbf{x}^{t+1} - (x^{t}+h\mathbf{v}+h^2\mathbf{M}^{-1}\mathbf{f}_{ext}) = h^2\mathbf{M}^{-1}\mathbf{f}_{int}(x_{n+1})
$$
上面的式子其实是一个最小化问题的驻点。
令$\mathbf{y} = x^{t}+h\mathbf{v}+h^2\mathbf{M}^{-1}\mathbf{f}_{ext}$
我们可以得到其变分形式：
$$\mathbf{x}_{n+1} = \min{\frac{1}{2h^2}| \mathbf{x} - \mathbf{y}|^{T}\mathbf{M}| \mathbf{x} - \mathbf{y}|+E(\mathbf{x})}$$
其中$E(\mathbf{x})$是我们的内力产生的弹性势能。
关于$E(\mathbf{x})$，就是PBD，XPBD，PD等文章的关键地方。如果$E(\mathbf{x})$的刚度无穷大，弹性势能必须满足，就是PBD的方法(下面的1)，如果最小化问题弹性势能和动能之间的平衡，$E(\mathbf{x})$刚度可控，具有特殊的形式那么就是XPBD（下面的2），PD方法（下面的3）。
- Gauss最小二乘问题
$$\min{\frac{1}{2h^2}| \mathbf{x} - \mathbf{y}|^{T}\mathbf{M}| \mathbf{x} - \mathbf{y}|} \\ s.t. C(\mathbf{x})= 0$$
这种方法的特点是约束必须满足，也就是硬约束。一般会使用拉格朗日乘子法来求解这个问题。 - 利用约束来定义弹性势能：$E(\mathbf{x}) = \frac{1}{2}C(\mathbf{x})^{T}\mathbf{K}C(\mathbf{x}) = \frac{1}{2}C(\mathbf{x})^{T}\mathbf{\alpha}^{-1}C(\mathbf{x})$， 其中$\mathbf{K}$为刚度矩阵，${\alpha}^{-1}$为服从矩阵。
$$\mathbf{x}_{n+1} = \min{\frac{1}{2h^2}| \mathbf{x} - \mathbf{y}|^{T}\mathbf{M}| \mathbf{x} - \mathbf{y}|+\frac{1}{2}C(\mathbf{x})^{T}\mathbf{\alpha}^{-1}C(\mathbf{x})}$$ - 利用特殊的能量表达式比如二次型表达式: $E(\mathbf{x}, \mathbf{p}) = d(\mathbf{x},\mathbf{p}) + \delta_{E}(\mathbf{p}) = \frac{w}{2}| \mathbf{A}\mathbf{x} - \mathbf{B}\mathbf{p}|^2_{F} + \delta_{E}(\mathbf{p})$。其中
$$\delta_{E}(\mathbf{p}) =
\left\{\begin{matrix}
0, & E(\mathbf{p}) =0,\\
+\infty, & otherwise.
\end{matrix}\right.$$
上面的方法2和方法3其实都是软约束，也就是约束在求解之后可能并没有满足。
XPBD求解方法
我们主要关注方法2中提到的最小化问题。在此之前我们应该了解两个知识点：
- 能量和力的关系
$$\mathbf{f}_{int} = -\nabla_{\mathbf{x}} E = - \nabla_{\mathbf{x}}(\frac{1}{2}C(\mathbf{x})^{T}\mathbf{\alpha}^{-1}C(\mathbf{x})) = -\nabla_{\mathbf{x}}C(\mathbf{x})\alpha^{-1}C(\mathbf{x})$$ - 我们假设$\lambda=\tilde{\alpha}^{-1}C(\mathbf{x}), \tilde{\alpha} = \frac{\alpha}{\Delta t^2}$,那么$\mathbf{f}_{int} = -\nabla_{\mathbf{x}}C(\mathbf{x})^T \lambda$。
所以我们最小化问题的解就是：
$$ \frac{\mathbf{M}}{h^2} (x-y) + \nabla_{\mathbf{x}}C(\mathbf{x})^T \lambda = 0$$
$$ \mathbf{c}(\mathbf{x}) - \alpha \lambda = 0$$
Strain Constraint
XPBD FEM
一个四面体的rest configuration下四个顶点的material coordinates$\mathbf{X}_0,\mathbf{X}_1,\mathbf{X}_2,\mathbf{X}_3$,变形后的world configuration下四个顶点的world coordinates为 $\mathbf{x}_0,\mathbf{x}_1,\mathbf{x}_2,\mathbf{x}_3$。我们可以得到矩阵
$$\mathbf{D}_m = [\mathbf{X}_1-\mathbf{X}_0, \mathbf{X}_2-\mathbf{X}_0,\mathbf{X}_3-\mathbf{X}_0],$$
令
$$
\mathbf{D}_m^{-1} = [\mathbf{c}_1,\mathbf{c}_2,\mathbf{c}_3],
$$
$$\mathbf{D}_w = [\mathbf{x}_1-\mathbf{x}_0,\mathbf{x}_2-\mathbf{x}_0,\mathbf{x}_3-\mathbf{x}_0]=[\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3].$$
那么变形梯度为：
$$\mathbf{F} = \mathbf{D}_w \mathbf{D}^{-1}_m = [\mathbf{f}_1,\mathbf{f}_2,\mathbf{f}_3],$$
我们可以得到Green Strain
$$\epsilon = \frac{1}{2} (\mathbf{F}^T \mathbf{F} - \mathbf{I}) = \frac{1}{2}(s_{ij} - \delta_{ij})= \frac{1}{2} [( \mathbf{D}_w \mathbf{c}_i ) \cdot (\mathbf{D}_w \mathbf{c}_j) - \delta_{ij} ].$$
Green strain 的矩阵形式如下：
$$\epsilon =
\begin{bmatrix}
\epsilon_x& 0.5\epsilon_{xy} & 0.5\epsilon_{xz}\\
0.5\epsilon_{xy}& \epsilon_y & 0.5\epsilon_{yz}\\
0.5\epsilon_{xz}& 0.5\epsilon_{yz}& \epsilon_z
\end{bmatrix},
$$
从上面可以看出$\epsilon$是对称的，也就是说有6个独立的分量，Voigt notation如下
$$ \mathbf{c} = [\epsilon_x,\epsilon_y,\epsilon_z,\epsilon_{xy},\epsilon_{yz},\epsilon_{zx}]^{T}= [s_{11}-1,s_{22}-1,s_{33}-1,s_{12},s_{13},s_{23}]^{T} \in \mathbb{R}^{6 \times 1}.$$
我们可以得到StVK模型的能量表达式为：
$$E(\mathbf{F}(\mathbf{x})) = V(\mu | \epsilon|^2_{F}+\frac{\lambda}{2}tr^2(\epsilon)) = \frac{1}{2}\mathbf{c}^T\mathbf{K}\mathbf{c}$$
其中 $V$是四面体rest configuration的体积，$\mu$ 和 $\lambda$ 是拉梅常数，$\mathbf{K} \in \mathbb{R}^{6 \times 6}$是刚度矩阵
$$\mathbf{K} = V
\begin{bmatrix}
\lambda + 2\mu& \lambda & \lambda & & & \\
\lambda& \lambda + 2\mu& \lambda& & & \\
\lambda& \lambda& \lambda + 2\mu& & & \\
& & & \mu& & \\
& & & & \mu& \\
& & & & & \mu
\end{bmatrix} = \mathbf{\alpha}^{-1}$$
在XPBD中，我们将$\mathbf{c}$作为约束进行求解。接下来需要进行两个步骤
计算Lagrangian乘子
$$
\Delta \lambda = \frac{-\mathbf{c} - \tilde{\alpha} \lambda}{\nabla \mathbf{c}^{T} \mathbf{M}^{-1}\nabla \mathbf{c} + \tilde{\alpha}}
$$
其中$\lambda \in \mathbb{R}^{6\times 1}$, $\tilde{\alpha} = \frac{\alpha}{h^2} \in \mathbb{R}^{6\times 6}, $h$是迭代步长, \mathbf{M} \in \mathbb{R}^{12\times12}$，它的对角线块矩阵是$m_0\mathbf{I}_{3\times3}，m_1\mathbf{I}_{3\times3}，m_2\mathbf{I}_{3\times3}，m_3\mathbf{I}_{3\times3}$。
要计算$\Delta \lambda$我们还需要计算约束$\mathbf{c}$分别关于$[ \mathbf{p}_0 , \mathbf{p}_1 , \mathbf{p}_2 , \mathbf{p}_3]$ 的梯度。
$$\nabla s_{ij} = \nabla [\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3] s_{ij}=\mathbf{f}_{i}\mathbf{c}_{j}^{T}+ \mathbf{f}_{j}\mathbf{c}_{i}^{T}$$
$$\nabla_{\mathbf{p}_0}s_{ij} = -\sum_{k=1}^{3}\nabla_{\mathbf{p}_k} s_{ij}$$
所以
$$\nabla \mathbf{c} =
\begin{bmatrix}
\frac{\partial\epsilon_x}{\partial\mathbf{p}_0} & \frac{\partial\epsilon_y}{\partial\mathbf{p}_0} & \frac{\partial\epsilon_z}{\partial\mathbf{p}_0} & \frac{\partial\epsilon_{xy}}{\partial\mathbf{p}_0} & \frac{\partial\epsilon_{yz}}{\partial\mathbf{p}_0} & \frac{\partial\epsilon_{zx}}{\partial\mathbf{p}_0} \\
\frac{\partial\epsilon_x}{\partial\mathbf{p}_1} & \frac{\partial\epsilon_y}{\partial\mathbf{p}_1} & \frac{\partial\epsilon_z}{\partial\mathbf{p}_1} & \frac{\partial\epsilon_{xy}}{\partial\mathbf{p}_1} & \frac{\partial\epsilon_{yz}}{\partial\mathbf{p}_1} & \frac{\partial\epsilon_{zx}}{\partial\mathbf{p}_1} \\
\frac{\partial\epsilon_x}{\partial\mathbf{p}_2} & \frac{\partial\epsilon_y}{\partial\mathbf{p}_2} & \frac{\partial\epsilon_z}{\partial\mathbf{p}_2} & \frac{\partial\epsilon_{xy}}{\partial\mathbf{p}_2} & \frac{\partial\epsilon_{yz}}{\partial\mathbf{p}_2} & \frac{\partial\epsilon_{zx}}{\partial\mathbf{p}_2} \\
\frac{\partial\epsilon_x}{\partial\mathbf{p}_3} & \frac{\partial\epsilon_y}{\partial\mathbf{p}_3} & \frac{\partial\epsilon_z}{\partial\mathbf{p}_3} & \frac{\partial\epsilon_{xy}}{\partial\mathbf{p}_3} & \frac{\partial\epsilon_{yz}}{\partial\mathbf{p}_3} & \frac{\partial\epsilon_{zx}}{\partial\mathbf{p}_3}
\end{bmatrix} =
\begin{bmatrix}
\nabla_{\mathbf{p}_0} s_{11} & \nabla_{\mathbf{p}_0}s_{22} & \nabla_{\mathbf{p}_0}s_{33} & \nabla_{\mathbf{p}_0} s_{12} & \nabla_{\mathbf{p}_0}s_{23} & \nabla_{\mathbf{p}_0}s_{13} \\
\nabla_{\mathbf{p}_1} s_{11} & \nabla_{\mathbf{p}_1}s_{22} & \nabla_{\mathbf{p}_1}s_{33} & \nabla_{\mathbf{p}_1} s_{12} & \nabla_{\mathbf{p}_1}s_{23} & \nabla_{\mathbf{p}_1}s_{13} \\
\nabla_{\mathbf{p}_2} s_{11} & \nabla_{\mathbf{p}_2}s_{22} & \nabla_{\mathbf{p}_2}s_{33} & \nabla_{\mathbf{p}_2} s_{12} & \nabla_{\mathbf{p}_2}s_{23} & \nabla_{\mathbf{p}_2}s_{13} \\
\nabla_{\mathbf{p}_3} s_{11} & \nabla_{\mathbf{p}_3}s_{22} & \nabla_{\mathbf{p}_3}s_{33} & \nabla_{\mathbf{p}_3} s_{12} & \nabla_{\mathbf{p}_3}s_{23} & \nabla_{\mathbf{p}_3}s_{13}
\end{bmatrix}
\in \mathbb{R}^{12 \times 6}
$$
计算修正位置
同理，根据公式$(17)$可以得到修正位置。
$$\lambda = \lambda + \Delta \lambda$$
$$
\Delta \mathbf{p}_{i} = \mathbf{m}_{i}^{-1} \nabla_{\mathbf{p}_{i}} \mathbf{c}_{i} \lambda
$$
其中$\nabla c_{i} \in \mathbb{R}^{3 \times 6}$
参考文献：
- 2014_SCA_Strain Based Dynamics
- 2016_SCA_XPBD: Position-Based Simulation of Compliant Constrained Dynamics
XPBD Rod
对零拉伸，可拉伸和扭转的Rod进行XPBD的建模：Rod离散成多个刚体，刚体之间利用零拉伸约束和拉伸扭转约束连接，具体的约束形式如下：
$$
\mathbf{C} (\mathbf{x}_1,q_1,\mathbf{x}_2,q_2) =
\begin{pmatrix}
\mathbf{R}(q_1)\mathbf{p}_{1}+\mathbf{x}_1 - \mathbf{R}(q_2)\mathbf{p}_{2}+\mathbf{x}_2 \\
\frac{2}{l_i}\Im [\bar{q}_1q_2-\bar{q}^{0}_{1}q^{0}_2]
\end{pmatrix} =
\begin{pmatrix}
\mathbf{C}_1 \\
\mathbf{C}_2
\end{pmatrix}
$$
有了约束之后，我们要对约束求梯度：
$$
\nabla \mathbf{C} =
\begin{pmatrix}
\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_1} & \frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_1} \\
\frac{\partial \mathbf{C}_1}{\partial q_1} & \frac{\partial \mathbf{C}_2}{\partial q_1} \\
\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_2} & \frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_2}\\
\frac{\partial \mathbf{C}_1}{\partial q_2} & \frac{\partial \mathbf{C}_2}{\partial q_2} \\
\end{pmatrix} \in \mathbb{R}^{12 \times 6}
$$
下面是具体的解析式：
$$
\begin{pmatrix}
\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_1} = \mathbf{1}_{3\times3} \\
\frac{\partial \mathbf{C}_1}{\partial q_1} = -[\mathbf{R}(q_1)\mathbf{p}_1]^{\times}\\
\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_2} = - \mathbf{1}_{3\times3} \\
\frac{\partial \mathbf{C}_1}{\partial q_2} = -[\mathbf{R}(q_2)\mathbf{p}_2]^{\times}\\
\end{pmatrix} \in \mathbb{R}^{12 \times 3},
\begin{pmatrix}
\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_1} = \mathbf{0}_{3\times3} \\
\frac{\partial \mathbf{C}_2}{\partial q_1} = \frac{\partial \Omega}{\partial q_1}\mathbf{G}(q_1)\\
\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_2} = - \mathbf{0}_{3\times3} \\
\frac{\partial \mathbf{C}_2}{\partial q_2} =\frac{\partial \Omega}{\partial q_2}\mathbf{G}(q_2)\\
\end{pmatrix} \in \mathbb{R}^{12 \times 3}
$$
接下来我们要求Lagrange 乘子的分母部分：$\nabla \mathbf{C}^{T}\mathbf{M} \nabla \mathbf{C}$，其中$\mathbf{M} \in \mathbb{R}^{12 \times 12}$ 具体表达式如下：
$$\mathbf{M} =
\begin{pmatrix}
m_1 \mathbf{1}_{3\times3} & && \\
& \mathbf{I}^{-1}_{1} &&& \\
&&m_2 \mathbf{1}_{3\times3} && \\
&&&\mathbf{I}^{-1}_{2}
\end{pmatrix}
$$
由此我们可以计算：
\begin{equation}
\nabla \mathbf{C}^{T}\mathbf{M} \nabla \mathbf{C} =
\begin{pmatrix}
\mathbf{A} & \mathbf{B} \\
\mathbf{B}^{T} & \mathbf{D}
\end{pmatrix}
\end{equation}
其中：
$$
\require{cancel}
\begin{aligned}
\mathbf{A} = &m_1\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_1} (\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_1})^{T} + \frac{\partial \mathbf{C}_1}{\partial q_1} \mathbf{I}^{-1}_{1} (\frac{\partial \mathbf{C}_1}{\partial q_1})^{T} +m_2\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_2} (\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_2})^{T} + \frac{\partial \mathbf{C}_1}{\partial q_2} \mathbf{I}^{-1}_{2} (\frac{\partial \mathbf{C}_1}{\partial q_2})^{T} \\
\mathbf{B} =& \cancel{ m_1\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_1} (\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_1})^{T}}+ \frac{\partial \mathbf{C}_1}{\partial q_1} \mathbf{I}^{-1}_{1} (\frac{\partial \mathbf{C}_2}{\partial q_1})^{T} + \cancel{m_2\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_2} (\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_2})^{T} }+ \frac{\partial \mathbf{C}_1}{\partial q_2} \mathbf{I}^{-1}_{2} (\frac{\partial \mathbf{C}_2}{\partial q_2})^{T} \\
\mathbf{B}^{T} =&\cancel{ m_1\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_1} (\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_1})^{T}}+ \frac{\partial \mathbf{C}_2}{\partial q_1} \mathbf{I}^{-1}_{1} (\frac{\partial \mathbf{C}_1}{\partial q_1})^{T} +\cancel{m_2\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_2} (\frac{\partial \mathbf{C}_1}{\partial \mathbf{x}_2})^{T}}+ \frac{\partial \mathbf{C}_2}{\partial q_2} \mathbf{I}^{-1}_{2} (\frac{\partial \mathbf{C}_1}{\partial q_2})^{T} \\
\mathbf{D} =& \cancel{m_1\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_1} (\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_1})^{T}}+ \frac{\partial \mathbf{C}_2}{\partial q_1} \mathbf{I}^{-1}_{1} (\frac{\partial \mathbf{C}_2}{\partial q_1})^{T} + \cancel{m_2\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_2} (\frac{\partial \mathbf{C}_2}{\partial \mathbf{x}_2})^{T}}+ \frac{\partial \mathbf{C}_2}{\partial q_2} \mathbf{I}^{-1}_{2} (\frac{\partial \mathbf{C}_2}{\partial q_2})^{T}
\end{aligned}
$$
将约束的具体解析式代入到$ \mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D}$中可以得到：
$$
\begin{aligned}
\mathbf{A} &= m_1 \mathbf{1}_{3\times3}+ [\mathbf{R}(q_1)\mathbf{p}_1]^{\times} \mathbf{I}^{-1}_{1} ([\mathbf{R}(q_1)\mathbf{p}_1]^{\times})^{T} +m_1 \mathbf{1}_{3\times3}+ [\mathbf{R}(q_2)\mathbf{p}_2]^{\times} \mathbf{I}^{-1}_{2} ( [\mathbf{R}(q_2)\mathbf{p}_2]^{\times} )^{T} \\
\mathbf{B} &= -[\mathbf{R}(q_1)\mathbf{p}_1]^{\times}\ \mathbf{I}^{-1}_{1} (\frac{\partial \Omega}{\partial q_1}\mathbf{G}(q_1))^{T} -[\mathbf{R}(q_2)\mathbf{p}_2] \mathbf{I}^{-1}_{2}(\frac{\partial \Omega}{\partial q_2}\mathbf{G}(q_2))^{T} \\
\mathbf{B}^{T} &= \mathbf{B}^{T}\\
\mathbf{D} &= \frac{\partial \Omega}{\partial q_1}\mathbf{G}(q_1) \mathbf{I}^{-1}_{1}(\frac{\partial \Omega}{\partial q_1}\mathbf{G}(q_1))^{T} + \frac{\partial \Omega}{\partial q_2}\mathbf{G}(q_2) \mathbf{I}^{-1}_{2}(\frac{\partial \Omega}{\partial q_2}\mathbf{G}(q_2))^{T}
\end{aligned}
$$
更新
更新拉格朗日乘子
$$\Delta \lambda = \frac{\mathbf{C} + \tilde{\alpha}\lambda}{ \nabla \mathbf{C}^{T}\mathbf{M} \nabla \mathbf{C} +\tilde{\alpha} } \in \mathbb{R}^{6\times1}$$
$$\lambda += \Delta \lambda$$
修正预测位置
$$(\Delta \mathbf{x}_i,\Delta q_i) =-w_i \nabla _{p_i,q_i}\mathbf{C} \lambda $$
参考文献
Direct Position-Based Solver for Stiff Rods
2D三角形网格基于XPBD的FEM （待完善）
对于一个二维的线性三角形单元，3个顶点分别是$\mathbf{p}_1,\mathbf{p}_2,\mathbf{p}_3$，对应的坐标$(x_1,y_1),(x_2,y_2),(x_3,y_3)$,对应的位移为$(\mu _1, v_1),(\mu _2, v_2), (\mu_3, v_3)$，令位移向量$\mathbf{q}= [\mu _1, v_1,\mu _2, v_2,\mu_3, v_3]^T$。
每个三角形单元就有6个节点位移即6个自由度。其对应的Strain $\varepsilon$为：
$$
\varepsilon =
\begin{Bmatrix}
\varepsilon_{x}\
\varepsilon_{y}\
\gamma_{xy}
\end{Bmatrix} =
\begin{Bmatrix}
\frac{\partial \mu}{\partial x}\\
\frac{\partial v}{\partial y}\\
\frac{\partial \mu}{\partial y} + \frac{\partial v}{\partial x}
\end{Bmatrix} = \mathbf{B}\mathbf{q}
$$
其中的$\mathbf{B} \in \mathbb{R}^{3\times 6}$ 是应变矩阵，具体的表达式如下：
$$
\mathbf{B} = \frac{1}{\det J}
\begin{bmatrix}
y_{23} & 0 & y_{31} & 0 & y_{12} & 0 \\
0 & x_{32} & 0 & x_{13} & 0 & x_{21}\\
x_{32} & y_{23} & x_{13} & y_{31} & x_{21} & y_{12}
\end{bmatrix}
$$
其中$y_{12} = y_1 - y_2$，其他同理。
Gradient of Straint Constraint
令 $\mathbf{B} = [\mathbf{B}_1, \mathbf{B}_2,\mathbf{B}_3]$,
则
$$
\varepsilon = [\mathbf{B}_1, \mathbf{B}_2,\mathbf{B}_3]
\begin{bmatrix}
\Delta \mathbf{p}_1\\
\Delta \mathbf{p}_2\\
\Delta \mathbf{p}_3
\end{bmatrix} =\mathbf{B}_1 \Delta \mathbf{p}_1+ \mathbf{B}_2 \Delta \mathbf{p}_2+ \mathbf{B}_3 \Delta \mathbf{p}_3
$$
令
$$\Delta \mathbf{p}_i = \mathbf{\tilde{p}}_i - \mathbf{p}_i$$
则
$$\varepsilon = \mathbf{B}_1 (\mathbf{\tilde{p}}_1 - \mathbf{p}_1))+ \mathbf{B}_2 (\mathbf{\tilde{p}}_2 - \mathbf{p}_2)+ \mathbf{B}_3 (\mathbf{\tilde{p}}_3 - \mathbf{p}_3)$$
对Strain $\varepsilon$关于 $\mathbf{p}_i$ 进行求导，可以得到对应的梯度
$$
\nabla \varepsilon_{\mathbf{p}_i} = \mathbf{B}^{T}_i
$$
最后得到总梯度为
$$
\nabla \varepsilon=[\mathbf{B}^{T}_1, \mathbf{B}^{T}_2,\mathbf{B}^{T}_3]^{T} \in \mathbb{R}^{6\times 3}
$$
代入XPBD中
计算Lagrangian乘子
$$
\Delta \lambda = \frac{-\varepsilon - \tilde{\alpha}_{tri} \lambda}{\nabla \varepsilon^{T} \mathbf{M}^{-1}\nabla \varepsilon + \tilde{\alpha}_{tri}}
$$
其中$\lambda \in \Re^{3\times 1}$, $\tilde{\alpha}_{tri} \in \Re^{3\times 3}, \mathbf{M} \in \Re^{6\times6}$。
计算修正位置
同理，根据公式$(17)$可以得到修正位置。
$$
\Delta \mathbf{p} = \mathbf{M}^{-1} \nabla \varepsilon \Delta \lambda
$$
具体推导步骤如下
形函数
我们以2D的三角形作为一个单元(Element)，3个顶点分别是$i,j,k$，对应的坐标$(x_i,y_i),(x_j,y_j),(x_k,y_k)$,对应的位移为$(\mu _i, v_i),(\mu _j, v_j), (\mu_k, v_k)$。
每个三角形单元就有6个节点位移即6个自由度。
我们假设每个位置$(x,y)$的位移$(\mu,v)$都是$(x,y)$的函数（可以想一想这是为什么？），进一步假设这个函数是一个多项式函数，随着项数和阶数的增加，这个多项式函数可以近似任意一个函数。在三角形单元中我们使用一阶的多项式来近似，即：
$$
\mu = \beta_1 + \beta_2 x+ \beta_3 y\
v = \beta_4 + \beta_5 x + \beta_6 y
$$