XPBD with Strain Energy

这一篇讲一讲怎么将Strain作为Constraint放到XPBD的框架里。

XPBD

在前面的博客当中我已经讲过XPBD的由来和概念了。现在来回顾一下推导的步骤:

  1. 由隐式欧拉法得到一个运动方程
  2. 得到运动方程的等价形式:最小值优化问题
  3. 将约束和能量等价,引入拉格朗日乘子
  4. 得到最小化问题的解(未简化)
  5. 简化处理得到问题的解

运动的连续时间积分

$$\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)。

  1. Gauss最小二乘问题
    $$\min{\frac{1}{2h^2}| \mathbf{x} - \mathbf{y}|^{T}\mathbf{M}| \mathbf{x} - \mathbf{y}|} \\ s.t. C(\mathbf{x})= 0$$
    这种方法的特点是约束必须满足,也就是硬约束。一般会使用拉格朗日乘子法来求解这个问题。
  2. 利用约束来定义弹性势能:$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})}$$
  3. 利用特殊的能量表达式比如二次型表达式: $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中提到的最小化问题。在此之前我们应该了解两个知识点:

  1. 能量和力的关系
    $$\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})$$
  2. 我们假设$\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}$

参考文献:

  1. 2014_SCA_Strain Based Dynamics
  2. 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
$$

具体推导步骤如下

形函数

Triangle Element
我们以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
$$

文章目录
  1. 1. XPBD
    1. 1.1. 运动的连续时间积分
    2. 1.2. 运动离散时间积分隐式欧拉法
    3. 1.3. 最小化问题
    4. 1.4. XPBD求解方法
  2. 2. Strain Constraint
    1. 2.1. XPBD FEM
      1. 2.1.1. 计算Lagrangian乘子
      2. 2.1.2. 计算修正位置
      3. 2.1.3. 参考文献:
  3. 3. XPBD Rod
    1. 3.1. 更新
      1. 3.1.1. 更新拉格朗日乘子
      2. 3.1.2. 修正预测位置
      3. 3.1.3. 参考文献
  4. 4. 2D三角形网格基于XPBD的FEM (待完善)
    1. 4.1. Gradient of Straint Constraint
    2. 4.2. 代入XPBD中
      1. 4.2.1. 计算Lagrangian乘子
      2. 4.2.2. 计算修正位置
    3. 4.3. 具体推导步骤如下
      1. 4.3.1. 形函数
,