我们已经知道力是势函数的梯度,但是在我们的离散化手段中,需要额外添加一个系数
f = − a ∇ Ψ ( x ) = − a ∇ ∂ Ψ ∂ x \mathbf{f}=-a\nabla \Psi(\mathbf{x})=-a\nabla\frac{\partial \Psi}{\partial \mathbf{x}} f = − a ∇Ψ ( x ) = − a ∇ ∂ x ∂ Ψ
a a a 是离散单元的度量,例如三角形的面积或空间四面体的体积,显然很小的三角形被压缩时不会和很大的三角形被压缩时具有同等的反馈力
我们还知道,不管是什么弹性势能函数 Ψ \Psi Ψ 都是根据变形梯度 F \mathbf{F} F 定义的,也就是以 F \mathbf{F} F 为参数,而 F \mathbf{F} F 又是根据离散元的静态坐标 (常量) x ‾ \mathbf{\overline{x}} x 和变形后的坐标 (变量) x \mathbf{x} x 定义的,所以计算力需要使用链式法则
∂ Ψ ∂ x = ∂ Ψ ∂ F ∂ F ∂ x \frac{\partial\Psi}{\partial\mathbf{x}}
=\frac{\partial\Psi}{\partial \mathbf{F}}\frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ x ∂ Ψ = ∂ F ∂ Ψ ∂ x ∂ F
我们先来计算 ∂ Ψ ∂ F \dfrac{\partial\Psi}{\partial\mathbf{F}} ∂ F ∂ Ψ ,在三维情况下,F \mathbf{F} F 是一个 3 × 3 3 \times 3 3 × 3 矩阵
F = [ f 00 f 01 f 02 f 10 f 11 f 12 f 20 f 21 f 22 ] = [ f 0 f 1 f 2 ] \mathbf{F} =
\begin{bmatrix}
f_{00} & f_{01} & f_{02} \\
f_{10} & f_{11} & f_{12} \\
f_{20} & f_{21} & f_{22}
\end{bmatrix}
=
\begin{bmatrix}
\mathbf{f}_0 & \mathbf{f}_1 & \mathbf{f}_2
\end{bmatrix} F = f 00 f 10 f 20 f 01 f 11 f 21 f 02 f 12 f 22 = [ f 0 f 1 f 2 ]
由于弹性势能函数 Ψ \Psi Ψ 是一个标量函数,而变形梯度 F \mathbf{F} F 是一个矩阵,所以偏导数 ∂ Ψ ∂ F \dfrac{\partial\Psi}{\partial\mathbf{F}} ∂ F ∂ Ψ 的结果也是一个矩阵
∂ Ψ ∂ F = [ ∂ Ψ ∂ f 00 ∂ Ψ ∂ f 01 ∂ Ψ ∂ f 02 ∂ Ψ ∂ f 10 ∂ Ψ ∂ f 11 ∂ Ψ ∂ f 12 ∂ Ψ ∂ f 20 ∂ Ψ ∂ f 21 ∂ Ψ ∂ f 22 ] ∈ R 3 × 3 \frac{\partial\Psi}{\partial\mathbf{F}} =
\left[\begin{array}{c c c}
\frac{\partial\Psi}{\partial f_{00}} & \frac{\partial\Psi}{\partial f_{01}} & \frac{\partial\Psi}{\partial f_{02}}
\\[2mm]
\frac{\partial\Psi}{\partial f_{10}} & \frac{\partial\Psi}{\partial f_{11}} & \frac{\partial\Psi}{\partial f_{12}}
\\[2mm]
\frac{\partial\Psi}{\partial f_{20}} & \frac{\partial\Psi}{\partial f_{21}} & \frac{\partial\Psi}{\partial f_{22}}
\end{array}\right]
\in \mathbb{R}^{3 \times 3} ∂ F ∂ Ψ = ∂ f 00 ∂ Ψ ∂ f 10 ∂ Ψ ∂ f 20 ∂ Ψ ∂ f 01 ∂ Ψ ∂ f 11 ∂ Ψ ∂ f 21 ∂ Ψ ∂ f 02 ∂ Ψ ∂ f 12 ∂ Ψ ∂ f 22 ∂ Ψ ∈ R 3 × 3
矩阵的每个元素 ∂ Ψ ∂ f i j \dfrac{\partial\Psi}{\partial f_{ij}} ∂ f ij ∂ Ψ 表示当 F \mathbf{F} F 的特定分量 f i j f_{ij} f ij 发生微小变化时,势能函数 Ψ \Psi Ψ 的变化率
计算 ∂ x ∂ F
接下来压力给到计算 ∂ F ∂ x \frac{\partial \mathbf{F}}{\partial \mathbf{x}} ∂ x ∂ F ,麻烦之处在于 F \mathbf{F} F 虽然是 x \mathbf{x} x 的函数,但是 F \mathbf{F} F 本身并不像 Ψ \Psi Ψ 一样是标量,而是一个矩阵,在三维情况下 F ∈ R 3 × 3 \mathbf{F} \in \mathbb{R}^{3\times 3} F ∈ R 3 × 3 ,我们要求它对 x \mathbf{x} x 的偏导数,但是 x \mathbf{x} x 本身是一个向量,在三维情况下,使用四面体来离散化,那么 x ∈ R 12 \mathbf{x} \in \mathbb{R}^{12} x ∈ R 12 (四面体有四个点,每个点在三维中有三个坐标),那么如果我们计算 F \mathbf{F} F 对 x \mathbf{x} x 的偏导数,则需要对 x \mathbf{x} x 的所有 12 个分量都计算一次,结果将是 12 个矩阵
例如假设我们对 x 0 x_0 x 0 计算偏导数:
∂ F ∂ x 0 = [ ∂ f 00 ∂ x 0 ∂ f 01 ∂ x 0 ∂ f 02 ∂ x 0 ∂ f 10 ∂ x 0 ∂ f 11 ∂ x 0 ∂ f 12 ∂ x 0 ∂ f 20 ∂ x 0 ∂ f 21 ∂ x 0 ∂ f 22 ∂ x 0 ] ∈ R 3 × 3 \frac{\partial \mathbf{F}}{\partial x_0} =
\left[\begin{array}{c c c}
\frac{\partial f_{00}}{\partial x_0} & \frac{\partial f_{01}}{\partial x_0} & \frac{\partial f_{02}}{\partial x_0}
\\[2mm]
\frac{\partial f_{10}}{\partial x_0} & \frac{\partial f_{11}}{\partial x_0} & \frac{\partial f_{12}}{\partial x_0}
\\[2mm]
\frac{\partial f_{20}}{\partial x_0} & \frac{\partial f_{21}}{\partial x_0} & \frac{\partial f_{22}}{\partial x_0}
\end{array}\right]
\in \mathbb{R}^{3 \times 3} ∂ x 0 ∂ F = ∂ x 0 ∂ f 00 ∂ x 0 ∂ f 10 ∂ x 0 ∂ f 20 ∂ x 0 ∂ f 01 ∂ x 0 ∂ f 11 ∂ x 0 ∂ f 21 ∂ x 0 ∂ f 02 ∂ x 0 ∂ f 12 ∂ x 0 ∂ f 22 ∈ R 3 × 3
这样的矩阵还有 12 个,并且它们必须要以整体的形式来考虑,为了方便的处理这种情况,我们建立 张量 Tensor 的概念,我们说 ∂ F ∂ x ∈ R 3 × 3 × 12 \frac{\partial \mathbf{F}}{\partial \mathbf{x}} \in \mathbb{R}^{3\times 3\times 12} ∂ x ∂ F ∈ R 3 × 3 × 12 是一个三阶张量 (可以认为需要几个维度来描述就是几阶张量),我们可以认为它是这样的:
∂ F ∂ x = [ [ ∂ F ∂ x 0 ] [ ∂ F ∂ x 1 ] [ ∂ F ∂ x 2 ] [ ∂ F ∂ x 3 ] ⋮ [ ∂ F ∂ x 11 ] ] \frac{\partial \mathbf{F}}{\partial \mathbf{x}} =
\left[\begin{array}{c c c}
\left[\begin{array}{c}
\frac{\partial \mathbf{F}}{\partial x_0}
\end{array}\right]
\\[2mm]
\left[\begin{array}{c}
\frac{\partial \mathbf{F}}{\partial x_1}
\end{array}\right]
\\[2mm]
\left[\begin{array}{c}
\frac{\partial \mathbf{F}}{\partial x_2}
\end{array}\right]
\\[2mm]
\left[\begin{array}{c}
\frac{\partial \mathbf{F}}{\partial x_3}
\end{array}\right]
\\[2mm]
\vdots
\\[2mm]
\left[\begin{array}{c}
\frac{\partial \mathbf{F}}{\partial x_{11}}
\end{array}\right]
\end{array}\right] ∂ x ∂ F = [ ∂ x 0 ∂ F ] [ ∂ x 1 ∂ F ] [ ∂ x 2 ∂ F ] [ ∂ x 3 ∂ F ] ⋮ [ ∂ x 11 ∂ F ]
其中每一个 ∂ F ∂ x ⋯ \dfrac{\partial \mathbf{F}}{\partial x_{\cdots}} ∂ x ⋯ ∂ F 都是一个 3 × 3 3 \times 3 3 × 3 矩阵
计算
∂ x ∂ Ψ = ∂ F ∂ Ψ ∂ x ∂ F
有了三阶张量 ∂ F ∂ x \frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ x ∂ F 以后,我们需要计算 ∂ Ψ ∂ F ∂ F ∂ x \frac{\partial\Psi}{\partial \mathbf{F}}\frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ F ∂ Ψ ∂ x ∂ F ,但是要如何计算矩阵和三阶张量的相乘? 我们首先要引入新的算子
双缩并 double-contraction 算子
双缩并 double-contraction 算子,使用符号 : : : ,它是向量内积 (点乘) 的推广,向量内积是
x ⋅ y = ⟨ x , y ⟩ = x T y = [ x 0 x 1 x 2 ] T [ y 0 y 1 y 2 ] = x 0 y 0 + x 1 y 1 + x 2 y 2 \mathbf{x} \cdot \mathbf{y} = \langle \mathbf{x}, \mathbf{y} \rangle
= \mathbf{x}^T \mathbf{y}
=
\begin{bmatrix}
x_0 \\
x_1 \\
x_2
\end{bmatrix}^T
\begin{bmatrix}
y_0 \\
y_1 \\
y_2
\end{bmatrix}
= x_0y_0 + x_1y_1 + x_2y_2 x ⋅ y = ⟨ x , y ⟩ = x T y = x 0 x 1 x 2 T y 0 y 1 y 2 = x 0 y 0 + x 1 y 1 + x 2 y 2
而矩阵的双缩并也就是矩阵对应的元素相乘之和,也称为 Frobenius 内积
A : B = ⟨ A , B ⟩ F = ∑ i j A i j B i j \mathbf{A} : \mathbf{B} = \langle \mathbf{A}, \mathbf{B} \rangle_{F}
= \sum_{ij} A_{ij}B_{ij} A : B = ⟨ A , B ⟩ F = ij ∑ A ij B ij
因为矩阵的 Frobenius 范数 就是矩阵自身与自身的 Frobenius 内积 的平方根
∥ A ∥ F = ⟨ A , A ⟩ F = ∑ i j A i j 2 \|\mathbf{A}\|_F = \langle \mathbf{A}, \mathbf{A} \rangle_{F}
= \sqrt{\sum_{ij} A_{ij}^2} ∥ A ∥ F = ⟨ A , A ⟩ F = ij ∑ A ij 2
双缩并算子要求两个矩阵的维度是相等的,否则将导致未定义行为
A : B = [ a 0 a 2 a 1 a 3 ] [ b 0 b 2 b 1 b 3 ] = a 0 b 0 + a 1 b 1 + a 2 b 2 + a 3 b 3 \mathbf{A} : \mathbf{B} =
\begin{bmatrix}
a_0 & a_2 \\
a_1 & a_3
\end{bmatrix}
\begin{bmatrix}
b_0 & b_2 \\
b_1 & b_3
\end{bmatrix}
= a_0b_0 + a_1b_1 + a_2b_2 + a_3b_3 A : B = [ a 0 a 1 a 2 a 3 ] [ b 0 b 1 b 2 b 3 ] = a 0 b 0 + a 1 b 1 + a 2 b 2 + a 3 b 3
定义张量的双缩并计算为
A : B = [ [ a 0 a 2 a 1 a 3 ] [ a 4 a 6 a 5 a 7 ] [ a 8 a 10 a 9 a 11 ] ] [ b 0 b 2 b 1 b 3 ] = [ a 0 b 0 + a 1 b 1 + a 2 b 2 + a 3 b 3 a 4 b 0 + a 5 b 1 + a 6 b 2 + a 7 b 3 a 8 b 0 + a 9 b 1 + a 10 b 2 + a 11 b 3 ] \mathbb{A} : \mathbf{B} =
\begin{bmatrix}
\begin{bmatrix}
a_0 & a_2 \\
a_1 & a_3 \\
\end{bmatrix}
\\
\\
\begin{bmatrix}
a_4 & a_6 \\
a_5 & a_7 \\
\end{bmatrix}
\\
\\
\begin{bmatrix}
a_8 & a_{10} \\
a_9 & a_{11} \\
\end{bmatrix}
\end{bmatrix}
\begin{bmatrix}
b_0 & b_2 \\
b_1 & b_3
\end{bmatrix}
=
\begin{bmatrix}
a_0b_0 + a_1b_1 + a_2b_2 + a_3b_3 \\
a_4b_0 + a_5b_1 + a_6b_2 + a_7b_3 \\
a_8b_0 + a_9b_1 + a_{10}b_2 + a_{11}b_3
\end{bmatrix} A : B = [ a 0 a 1 a 2 a 3 ] [ a 4 a 5 a 6 a 7 ] [ a 8 a 9 a 10 a 11 ] [ b 0 b 1 b 2 b 3 ] = a 0 b 0 + a 1 b 1 + a 2 b 2 + a 3 b 3 a 4 b 0 + a 5 b 1 + a 6 b 2 + a 7 b 3 a 8 b 0 + a 9 b 1 + a 10 b 2 + a 11 b 3
所以 ∂ Ψ ∂ F ∂ F ∂ x \frac{\partial\Psi}{\partial \mathbf{F}}\frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ F ∂ Ψ ∂ x ∂ F 可以看作一种双缩并计算
∂ F ∂ x : ∂ Ψ ∂ F \frac{\partial\mathbf{F}}{\partial \mathbf{x}} : \frac{\partial\Psi}{\partial \mathbf{F}} ∂ x ∂ F : ∂ F ∂ Ψ
展平张量 Flattened Tensors
另一种计算矩阵和张量相乘的观点是进行展平
矩阵的展平是
vec ( A ) = vec ( [ a 0 a 2 a 1 a 3 ] ) = [ a 0 a 1 a 2 a 3 ] \text{vec}(\mathbf{A})
= \text{vec}
\left(
\begin{bmatrix}
a_0 & a_2 \\
a_1 & a_3 \\
\end{bmatrix}
\right)
=
\begin{bmatrix}
a_0 \\
a_1 \\
a_2 \\
a_3
\end{bmatrix} vec ( A ) = vec ( [ a 0 a 1 a 2 a 3 ] ) = a 0 a 1 a 2 a 3
张量的展平是
vec ( A ) = vec [ [ A ] [ B ] [ C ] ] = [ vec ( A ) vec ( B ) vec ( C ) ] \text{vec}(\mathbb{A}) = \text{vec}
\begin{bmatrix}
\begin{bmatrix}
\mathbf{A}
\end{bmatrix}
\\
\\
\begin{bmatrix}
\mathbf{B}
\end{bmatrix}
\\
\\
\begin{bmatrix}
\mathbf{C}
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\text{vec}(\mathbf{A}) & \text{vec}(\mathbf{B}) & \text{vec}(\mathbf{C})
\end{bmatrix} vec ( A ) = vec [ A ] [ B ] [ C ] = [ vec ( A ) vec ( B ) vec ( C ) ]
例如
vec ( A ) = vec [ [ a 0 a 2 a 1 a 3 ] [ a 4 a 6 a 5 a 7 ] [ a 8 a 10 a 9 a 11 ] ] = [ vec ( [ a 0 a 2 a 1 a 3 ] ) vec ( [ a 4 a 6 a 5 a 7 ] ) vec ( [ a 8 a 10 a 9 a 11 ] ) ] = [ a 0 a 4 a 8 a 1 a 5 a 9 a 2 a 6 a 10 a 3 a 7 a 11 ] \text{vec}(\mathbb{A}) = \text{vec}
\begin{bmatrix}
\begin{bmatrix}
a_0 & a_2 \\
a_1 & a_3 \\
\end{bmatrix}
\\
\\
\begin{bmatrix}
a_4 & a_6 \\
a_5 & a_7 \\
\end{bmatrix}
\\
\\
\begin{bmatrix}
a_8 & a_{10} \\
a_9 & a_{11} \\
\end{bmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\text{vec}\left(
\begin{bmatrix}
a_0 & a_2 \\
a_1 & a_3 \\
\end{bmatrix}
\right) &
\text{vec}\left(
\begin{bmatrix}
a_4 & a_6 \\
a_5 & a_7 \\
\end{bmatrix}
\right) &
\text{vec}\left(
\begin{bmatrix}
a_8 & a_{10} \\
a_9 & a_{11} \\
\end{bmatrix}
\right)
\end{bmatrix}
=
\begin{bmatrix}
a_0 & a_4 & a_8 \\
a_1 & a_5 & a_9 \\
a_2 & a_6 & a_{10} \\
a_3 & a_7 & a_{11}
\end{bmatrix} vec ( A ) = vec [ a 0 a 1 a 2 a 3 ] [ a 4 a 5 a 6 a 7 ] [ a 8 a 9 a 10 a 11 ] = [ vec ( [ a 0 a 1 a 2 a 3 ] ) vec ( [ a 4 a 5 a 6 a 7 ] ) vec ( [ a 8 a 9 a 10 a 11 ] ) ] = a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 a 10 a 11
这样一来,我们可以以展平形式得到和 双缩并 一致的结果
vec A T vec B = vec ( [ [ a 0 a 2 a 1 a 3 ] [ a 4 a 6 a 5 a 7 ] [ a 8 a 10 a 9 a 11 ] ] ) T vec [ b 0 b 2 b 1 b 3 ] = [ a 0 a 4 a 8 a 1 a 5 a 9 a 2 a 6 a 10 a 3 a 7 a 11 ] T [ b 0 b 1 b 2 b 3 ] = [ a 0 b 0 + a 1 b 1 + a 2 b 2 + a 3 b 3 a 4 b 0 + a 5 b 1 + a 6 b 2 + a 7 b 3 a 8 b 0 + a 9 b 1 + a 10 b 2 + a 11 b 3 ] \begin{align*}
\text{vec}\mathbb{A}^T \text{vec}\mathbf{B}
&=
\text{vec}
\left(
\begin{bmatrix}
\begin{bmatrix}
a_0 & a_2 \\
a_1 & a_3 \\
\end{bmatrix}
\\
\\
\begin{bmatrix}
a_4 & a_6 \\
a_5 & a_7 \\
\end{bmatrix}
\\
\\
\begin{bmatrix}
a_8 & a_{10} \\
a_9 & a_{11} \\
\end{bmatrix}
\end{bmatrix}
\right)^T
\text{vec}
\begin{bmatrix}
b_0 & b_2 \\
b_1 & b_3
\end{bmatrix}
\\
&=
\begin{bmatrix}
a_0 & a_4 & a_8 \\
a_1 & a_5 & a_9 \\
a_2 & a_6 & a_{10} \\
a_3 & a_7 & a_{11}
\end{bmatrix}^T
\begin{bmatrix}
b_0 \\
b_1 \\
b_2 \\
b_3
\end{bmatrix}\\
&=
\begin{bmatrix}
a_0b_0 + a_1b_1 + a_2b_2 + a_3b_3 \\
a_4b_0 + a_5b_1 + a_6b_2 + a_7b_3 \\
a_8b_0 + a_9b_1 + a_{10}b_2 + a_{11}b_3
\end{bmatrix}
\end{align*} vec A T vec B = vec [ a 0 a 1 a 2 a 3 ] [ a 4 a 5 a 6 a 7 ] [ a 8 a 9 a 10 a 11 ] T vec [ b 0 b 1 b 2 b 3 ] = a 0 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9 a 10 a 11 T b 0 b 1 b 2 b 3 = a 0 b 0 + a 1 b 1 + a 2 b 2 + a 3 b 3 a 4 b 0 + a 5 b 1 + a 6 b 2 + a 7 b 3 a 8 b 0 + a 9 b 1 + a 10 b 2 + a 11 b 3
最终计算力
根据前文我们已经指定力的公式是
f = − a ∂ Ψ ∂ x = ∂ Ψ ∂ F ∂ F ∂ x = ∂ F ∂ x : ∂ Ψ ∂ F \mathbf{f} = -a\frac{\partial\Psi}{\partial\mathbf{x}}
= \frac{\partial\Psi}{\partial \mathbf{F}}\frac{\partial\mathbf{F}}{\partial \mathbf{x}}
= \frac{\partial\mathbf{F}}{\partial \mathbf{x}} : \frac{\partial\Psi}{\partial \mathbf{F}} f = − a ∂ x ∂ Ψ = ∂ F ∂ Ψ ∂ x ∂ F = ∂ x ∂ F : ∂ F ∂ Ψ
可以看到,三阶张量项 ∂ F ∂ x \frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ x ∂ F 与能量函数 Ψ \Psi Ψ 无关,所以我们只需要实现它一次,然后对于不同的弹性势能函数 Ψ \Psi Ψ 去实现不同的 ∂ Ψ ∂ F \frac{\partial\Psi}{\partial \mathbf{F}} ∂ F ∂ Ψ 即可
这里指的三阶张量不变性是对于相同的离散元而言的,如果从四面体切换到六面体,那么变形梯度 F \mathbf{F} F 的结构显然会发生变化,那么 ∂ F ∂ x \frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ x ∂ F 也必然发生变化
计算 ∂ F ∂ x \frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ x ∂ F 张量
我们已经知道变形梯度是
x = [ x 0 x 1 x 2 x 3 ] x ‾ = [ x ‾ 0 x ‾ 1 x ‾ 2 x ‾ 3 ] F = D s D m − 1 = [ ∣ ∣ ∣ x 1 − x 0 x 2 − x 0 x 3 − x 0 ∣ ∣ ∣ ] [ ∣ ∣ ∣ x ‾ 1 − x ‾ 0 x ‾ 2 − x ‾ 0 x ‾ 3 − x ‾ 0 ∣ ∣ ∣ ] − 1 \begin{align*}
\mathbf{x} &=
\begin{bmatrix}
\mathbf{x}_0 \\
\mathbf{x}_1 \\
\mathbf{x}_2 \\
\mathbf{x}_3
\end{bmatrix}
\quad
\mathbf{\overline{x}} =
\begin{bmatrix}
\mathbf{\overline{x}}_0 \\
\mathbf{\overline{x}}_1 \\
\mathbf{\overline{x}}_2 \\
\mathbf{\overline{x}}_3
\end{bmatrix}
\\
\mathbf{F} &= \mathbf{D}_s \mathbf{D}_m^{-1}
\\
&=
\begin{bmatrix}
\ | &| &| \\
\mathbf{x_1} - \mathbf{x_0} & \mathbf{x_2} - \mathbf{x_0} & \mathbf{x_3} - \mathbf{x_0}\\
\ | &| &| \\
\end{bmatrix}
\begin{bmatrix}
\ | &| &| \\
\mathbf{\overline{x}_1} - \mathbf{\overline{x}_0} & \mathbf{\overline{x}_2} - \mathbf{\overline{x}_0} & \mathbf{\overline{x}_3} - \mathbf{\overline{x}_0}\\
\ | &| &| \\
\end{bmatrix}^{-1}
\end{align*} x F = x 0 x 1 x 2 x 3 x = x 0 x 1 x 2 x 3 = D s D m − 1 = ∣ x 1 − x 0 ∣ ∣ x 2 − x 0 ∣ ∣ x 3 − x 0 ∣ ∣ x 1 − x 0 ∣ ∣ x 2 − x 0 ∣ ∣ x 3 − x 0 ∣ − 1
所幸我们要求的是关于形变坐标 x \mathbf{x} x 的偏导,而非静态坐标 x ‾ \mathbf{\overline{x}} x ,那会需要处理矩阵逆的导数 (很麻烦),我们最终会得到这样的式子
∂ F ∂ x i = ∂ D s ∂ x i D m − 1 \frac{\partial\mathbf{F}}{\partial x_i}
=
\frac{\partial\mathbf{D}_s}{\partial x_i}\mathbf{D}_m^{-1} ∂ x i ∂ F = ∂ x i ∂ D s D m − 1
其中 D m − 1 \mathbf{D}_m^{-1} D m − 1 是以静态坐标 x ‾ \mathbf{\overline{x}} x 为输入参数的,所以只需要在一开始计算一次逆矩阵,然后如上所述的 ∂ D s ∂ x i \frac{\partial\mathbf{D}_s}{\partial x_i} ∂ x i ∂ D s 在三维时使用 四面体 tetrahedron 离散元时会产生 12 个矩阵组成的张量
我们已知形变坐标向量是
x = [ x 0 x 1 x 2 x 3 ] = [ x 0 x 1 ⋮ x 11 ] \mathbf{x} =
\begin{bmatrix}
\mathbf{x}_0 \\
\mathbf{x}_1 \\
\mathbf{x}_2 \\
\mathbf{x}_3
\end{bmatrix}
=
\begin{bmatrix}
x_0\\x_1\\ \vdots \\ x_{11}
\end{bmatrix} x = x 0 x 1 x 2 x 3 = x 0 x 1 ⋮ x 11
所以张量 ∂ D s ∂ x \frac{\partial\mathbf{D}_s}{\partial \mathbf{x}} ∂ x ∂ D s 是
∂ F ∂ x = [ [ ∂ D s ∂ x 0 ] [ ∂ D s ∂ x 1 ] [ ∂ D s ∂ x 2 ] [ ∂ D s ∂ x 3 ] ⋮ [ ∂ D s ∂ x 11 ] ] \frac{\partial \mathbf{F}}{\partial \mathbf{x}} =
\left[\begin{array}{c c c}
\left[\begin{array}{c}
\frac{\partial \mathbf{D}_s}{\partial x_0}
\end{array}\right]
\\[2mm]
\left[\begin{array}{c}
\frac{\partial \mathbf{D}_s}{\partial x_1}
\end{array}\right]
\\[2mm]
\left[\begin{array}{c}
\frac{\partial \mathbf{D}_s}{\partial x_2}
\end{array}\right]
\\[2mm]
\left[\begin{array}{c}
\frac{\partial \mathbf{D}_s}{\partial x_3}
\end{array}\right]
\\[2mm]
\vdots
\\[2mm]
\left[\begin{array}{c}
\frac{\partial \mathbf{D}_s}{\partial x_{11}}
\end{array}\right]
\end{array}\right] ∂ x ∂ F = [ ∂ x 0 ∂ D s ] [ ∂ x 1 ∂ D s ] [ ∂ x 2 ∂ D s ] [ ∂ x 3 ∂ D s ] ⋮ [ ∂ x 11 ∂ D s ]
我们已知四面体各顶点形变坐标
x 0 = [ x 0 x 1 x 2 ] x 1 = [ x 3 x 4 x 5 ] x 2 = [ x 6 x 7 x 8 ] x 3 = [ x 9 x 10 x 11 ] \mathbf{x_0} =
\begin{bmatrix}
x_0 \\
x_1 \\
x_2
\end{bmatrix}
\quad
\mathbf{x_1} =
\begin{bmatrix}
x_3 \\
x_4 \\
x_5
\end{bmatrix}
\quad
\mathbf{x_2} =
\begin{bmatrix}
x_6 \\
x_7 \\
x_8
\end{bmatrix}
\quad
\mathbf{x_3} =
\begin{bmatrix}
x_9 \\
x_{10} \\
x_{11}
\end{bmatrix} x 0 = x 0 x 1 x 2 x 1 = x 3 x 4 x 5 x 2 = x 6 x 7 x 8 x 3 = x 9 x 10 x 11
所以 D s \mathbf{D}_s D s 是
D s = [ ∣ ∣ ∣ x 1 − x 0 x 2 − x 0 x 3 − x 0 ∣ ∣ ∣ ] = [ x 3 − x 0 x 6 − x 0 x 9 − x 0 x 4 − x 1 x 7 − x 1 x 10 − x 1 x 5 − x 2 x 8 − x 2 x 11 − x 2 ] \begin{align*}
\mathbf{D}_s &=
\begin{bmatrix}
\ | &| &| \\
\mathbf{x_1} - \mathbf{x_0} & \mathbf{x_2} - \mathbf{x_0} & \mathbf{x_3} - \mathbf{x_0}\\
\ | &| &| \\
\end{bmatrix}
\\
&=
\begin{bmatrix}
x_3 - x_0 & x_6 - x_0 & x_9 - x_0 \\
x_4 - x_1 & x_7 - x_1 & x_{10} - x_1 \\
x_5 - x_2 & x_8 - x_2 & x_{11} - x_2
\end{bmatrix}
\end{align*} D s = ∣ x 1 − x 0 ∣ ∣ x 2 − x 0 ∣ ∣ x 3 − x 0 ∣ = x 3 − x 0 x 4 − x 1 x 5 − x 2 x 6 − x 0 x 7 − x 1 x 8 − x 2 x 9 − x 0 x 10 − x 1 x 11 − x 2
由此,我们可以计算张量 ∂ D s ∂ x i \frac{\partial\mathbf{D}_s}{\partial x_i} ∂ x i ∂ D s 各分量
∂ D s ∂ x 0 = [ − 1 − 1 − 1 0 0 0 0 0 0 ] ∂ D s ∂ x 1 = [ 0 0 0 − 1 − 1 − 1 0 0 0 ] ∂ D s ∂ x 2 = [ 0 0 0 0 0 0 − 1 − 1 − 1 ] ∂ D s ∂ x 3 = [ 1 0 0 0 0 0 0 0 0 ] ∂ D s ∂ x 4 = [ 0 0 0 1 0 0 0 0 0 ] ∂ D s ∂ x 5 = [ 0 0 0 0 0 0 1 0 0 ] ∂ D s ∂ x 6 = [ 0 1 0 0 0 0 0 0 0 ] ∂ D s ∂ x 7 = [ 0 0 0 0 1 0 0 0 0 ] ∂ D s ∂ x 8 = [ 0 0 0 0 0 0 0 1 0 ] ∂ D s ∂ x 9 = [ 0 0 1 0 0 0 0 0 0 ] ∂ D s ∂ x 10 = [ 0 0 0 0 0 1 0 0 0 ] ∂ D s ∂ x 11 = [ 0 0 0 0 0 0 0 0 1 ] \begin{array}{l l l}
\frac{\partial\mathbf{D}_s}{\partial x_0} =
\begin{bmatrix}
-1 & -1 & -1 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{D}_s}{\partial x_1} =
\begin{bmatrix}
0 & 0 & 0 \\
-1 & -1 & -1 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{D}_s}{\partial x_2} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
-1 & -1 & -1
\end{bmatrix}
\\[8mm]
\frac{\partial\mathbf{D}_s}{\partial x_3} =
\begin{bmatrix}
1 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{D}_s}{\partial x_4} =
\begin{bmatrix}
0 & 0 & 0 \\
1 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{D}_s}{\partial x_5} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
1 & 0 & 0
\end{bmatrix}
\\[8mm]
\frac{\partial\mathbf{D}_s}{\partial x_6} =
\begin{bmatrix}
0 & 1 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{D}_s}{\partial x_7} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{D}_s}{\partial x_8} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 1 & 0
\end{bmatrix}
\\[8mm]
\frac{\partial\mathbf{D}_s}{\partial x_9} =
\begin{bmatrix}
0 & 0 & 1 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{D}_s}{\partial x_{10}} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 1 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{D}_s}{\partial x_{11}} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 1
\end{bmatrix}
\end{array} ∂ x 0 ∂ D s = − 1 0 0 − 1 0 0 − 1 0 0 ∂ x 3 ∂ D s = 1 0 0 0 0 0 0 0 0 ∂ x 6 ∂ D s = 0 0 0 1 0 0 0 0 0 ∂ x 9 ∂ D s = 0 0 0 0 0 0 1 0 0 ∂ x 1 ∂ D s = 0 − 1 0 0 − 1 0 0 − 1 0 ∂ x 4 ∂ D s = 0 1 0 0 0 0 0 0 0 ∂ x 7 ∂ D s = 0 0 0 0 1 0 0 0 0 ∂ x 10 ∂ D s = 0 0 0 0 0 0 0 1 0 ∂ x 2 ∂ D s = 0 0 − 1 0 0 − 1 0 0 − 1 ∂ x 5 ∂ D s = 0 0 1 0 0 0 0 0 0 ∂ x 8 ∂ D s = 0 0 0 0 0 1 0 0 0 ∂ x 11 ∂ D s = 0 0 0 0 0 0 0 0 1
有了 ∂ D s ∂ x i \frac{\partial\mathbf{D}_s}{\partial x_i} ∂ x i ∂ D s 之后,我们要计算
∂ F ∂ x i = ∂ D s ∂ x i D m − 1 \frac{\partial\mathbf{F}}{\partial x_i}
=
\frac{\partial\mathbf{D}_s}{\partial x_i}\mathbf{D}_m^{-1} ∂ x i ∂ F = ∂ x i ∂ D s D m − 1
注意这里不是双缩并计算,因为张量 ∂ D s ∂ x i \frac{\partial\mathbf{D}_s}{\partial x_i} ∂ x i ∂ D s 不是我们想要的最终张量,我们需要对其中的每个矩阵乘以 D m − 1 \mathbf{D}_m^{-1} D m − 1 ,也就是说目前还处于生成最终张量 ∂ F ∂ x \frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ x ∂ F 的中间步骤
为了计算这些矩阵乘法,我们对 D m − 1 \mathbf{D}_m^{-1} D m − 1 定义一些有用的辅助记号
D m − 1 = [ d 00 d 01 d 02 d 10 d 11 d 12 d 20 d 21 d 22 ] r 0 = [ d 00 d 01 d 02 ] r 1 = [ d 10 d 11 d 12 ] r 2 = [ d 20 d 21 d 22 ] c 0 = [ d 00 d 10 d 20 ] T c 1 = [ d 01 d 11 d 21 ] T c 2 = [ d 02 d 12 d 22 ] T s 0 = d 00 + d 10 + d 20 s 1 = d 01 + d 11 + d 21 s 2 = d 02 + d 12 + d 22 D m − 1 = [ – r 0 – – r 1 – – r 2 – ] = [ ∣ ∣ ∣ c 0 c 1 c 2 ∣ ∣ ∣ ] \begin{align*}
\mathbf{D}_m^{-1} &=
\begin{bmatrix}
d_{00} & d_{01} & d_{02} \\
d_{10} & d_{11} & d_{12} \\
d_{20} & d_{21} & d_{22}
\end{bmatrix}
\\
\mathbf{r}_0 &= \begin{bmatrix} d_{00} & d_{01} & d_{02} \end{bmatrix}
\\
\mathbf{r}_1 &= \begin{bmatrix} d_{10} & d_{11} & d_{12} \end{bmatrix}
\\
\mathbf{r}_2 &= \begin{bmatrix} d_{20} & d_{21} & d_{22} \end{bmatrix}
\\ \\
\mathbf{c}_0 &= \begin{bmatrix} d_{00} & d_{10} & d_{20} \end{bmatrix}^T
\\
\mathbf{c}_1 &= \begin{bmatrix} d_{01} & d_{11} & d_{21} \end{bmatrix}^T
\\
\mathbf{c}_2 &= \begin{bmatrix} d_{02} & d_{12} & d_{22} \end{bmatrix}^T
\\ \\
s_0 &= d_{00} + d_{10} + d_{20}
\\
s_1 &= d_{01} + d_{11} + d_{21}
\\
s_2 &= d_{02} + d_{12} + d_{22}
\\ \\
\mathbf{D}_m^{-1} &=
\begin{bmatrix}
\text{--} \ \ \mathbf{r}_0 \ \text{--} \\
\text{--} \ \ \mathbf{r}_1 \ \text{--} \\
\text{--} \ \ \mathbf{r}_2 \ \text{--}
\end{bmatrix}
=
\begin{bmatrix}
| & | & | \\
\mathbf{c}_0 & \mathbf{c}_1 & \mathbf{c}_2 \\
| & | & | \\
\end{bmatrix}
\end{align*} D m − 1 r 0 r 1 r 2 c 0 c 1 c 2 s 0 s 1 s 2 D m − 1 = d 00 d 10 d 20 d 01 d 11 d 21 d 02 d 12 d 22 = [ d 00 d 01 d 02 ] = [ d 10 d 11 d 12 ] = [ d 20 d 21 d 22 ] = [ d 00 d 10 d 20 ] T = [ d 01 d 11 d 21 ] T = [ d 02 d 12 d 22 ] T = d 00 + d 10 + d 20 = d 01 + d 11 + d 21 = d 02 + d 12 + d 22 = – r 0 – – r 1 – – r 2 – = ∣ c 0 ∣ ∣ c 1 ∣ ∣ c 2 ∣
根据这些记号,我们可以计算出每个矩阵乘法结果
∂ F ∂ x 0 = [ − s 0 − s 1 − s 2 0 0 0 0 0 0 ] ∂ F ∂ x 1 = [ 0 0 0 − s 0 − s 1 − s 2 0 0 0 ] ∂ F ∂ x 2 = [ 0 0 0 0 0 0 − s 0 − s 1 − s 2 ] ∂ F ∂ x 3 = [ r 0 0 0 0 0 0 0 ] ∂ F ∂ x 4 = [ 0 0 0 r 0 0 0 0 ] ∂ F ∂ x 5 = [ 0 0 0 0 0 0 r 0 ] ∂ F ∂ x 6 = [ r 1 0 0 0 0 0 0 ] ∂ F ∂ x 7 = [ 0 0 0 r 1 0 0 0 ] ∂ F ∂ x 8 = [ 0 0 0 0 0 0 r 1 ] ∂ F ∂ x 9 = [ r 2 0 0 0 0 0 0 ] ∂ F ∂ x 10 = [ 0 0 0 r 2 0 0 0 ] ∂ F ∂ x 11 = [ 0 0 0 0 0 0 r 2 ] \begin{array}{l l l}
\frac{\partial\mathbf{F}}{\partial x_0} =
\begin{bmatrix}
-s_0 & -s_1 & -s_2 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{F}}{\partial x_1} =
\begin{bmatrix}
0 & 0 & 0 \\
-s_0 & -s_1 & -s_2 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{F}}{\partial x_2} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
-s_0 & -s_1 & -s_2
\end{bmatrix}
\\[8mm]
\frac{\partial\mathbf{F}}{\partial x_3} =
\begin{bmatrix}
&\mathbf{r}_0& \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{F}}{\partial x_4} =
\begin{bmatrix}
0 & 0 & 0 \\
&\mathbf{r}_0& \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{F}}{\partial x_5} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
&\mathbf{r}_0&
\end{bmatrix}
\\[8mm]
\frac{\partial\mathbf{F}}{\partial x_6} =
\begin{bmatrix}
&\mathbf{r}_1& \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{F}}{\partial x_7} =
\begin{bmatrix}
0 & 0 & 0 \\
&\mathbf{r}_1& \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{F}}{\partial x_8} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
&\mathbf{r}_1&
\end{bmatrix}
\\[8mm]
\frac{\partial\mathbf{F}}{\partial x_9} =
\begin{bmatrix}
&\mathbf{r}_2& \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{F}}{\partial x_{10}} =
\begin{bmatrix}
0 & 0 & 0 \\
&\mathbf{r}_2& \\
0 & 0 & 0
\end{bmatrix}
&
\frac{\partial\mathbf{F}}{\partial x_{11}} =
\begin{bmatrix}
0 & 0 & 0 \\
0 & 0 & 0 \\
&\mathbf{r}_2&
\end{bmatrix}
\end{array} ∂ x 0 ∂ F = − s 0 0 0 − s 1 0 0 − s 2 0 0 ∂ x 3 ∂ F = 0 0 r 0 0 0 0 0 ∂ x 6 ∂ F = 0 0 r 1 0 0 0 0 ∂ x 9 ∂ F = 0 0 r 2 0 0 0 0 ∂ x 1 ∂ F = 0 − s 0 0 0 − s 1 0 0 − s 2 0 ∂ x 4 ∂ F = 0 0 0 r 0 0 0 0 ∂ x 7 ∂ F = 0 0 0 r 1 0 0 0 ∂ x 10 ∂ F = 0 0 0 r 2 0 0 0 ∂ x 2 ∂ F = 0 0 − s 0 0 0 − s 1 0 0 − s 2 ∂ x 5 ∂ F = 0 0 0 0 r 0 0 0 ∂ x 8 ∂ F = 0 0 0 0 r 1 0 0 ∂ x 11 ∂ F = 0 0 0 0 r 2 0 0
最终结果都来自 D m − 1 \mathbf{D}_m^{-1} D m − 1 ,这 12 个矩阵最终组成了张量 ∂ F ∂ x \frac{\partial\mathbf{F}}{\partial \mathbf{x}} ∂ x ∂ F
如果此时我们还对某个弹性势能函数计算了 ∂ Ψ ∂ F \frac{\partial\Psi}{\partial \mathbf{F}} ∂ F ∂ Ψ
那么我们得以最终计算双缩并
∂ F ∂ x : ∂ Ψ ∂ F \frac{\partial\mathbf{F}}{\partial \mathbf{x}} : \frac{\partial\Psi}{\partial \mathbf{F}} ∂ x ∂ F : ∂ F ∂ Ψ
或者我们也可以用展平的方式做等价的计算
∂ Ψ ∂ x = ∂ Ψ ∂ F ∂ F ∂ x = ∂ F ∂ x : ∂ Ψ ∂ F = vec ( ∂ F ∂ x ) T vec ( ∂ Ψ ∂ F ) \frac{\partial\Psi}{\partial\mathbf{x}}
= \frac{\partial\Psi}{\partial \mathbf{F}}\frac{\partial\mathbf{F}}{\partial \mathbf{x}}
= \frac{\partial\mathbf{F}}{\partial \mathbf{x}} : \frac{\partial\Psi}{\partial \mathbf{F}}
= \text{vec}\left(\frac{\partial\mathbf{F}}{\partial \mathbf{x}} \right)^T
\text{vec}\left(\frac{\partial\Psi}{\partial \mathbf{F}} \right) ∂ x ∂ Ψ = ∂ F ∂ Ψ ∂ x ∂ F = ∂ x ∂ F : ∂ F ∂ Ψ = vec ( ∂ x ∂ F ) T vec ( ∂ F ∂ Ψ )