FEM中的梯度计算公式证明与matlab验证

证明

首先，已知

\[ \vec {\mathbf{x}} = \begin{bmatrix} \vec x_1 & \vec x_2 & \vec x_3 \end{bmatrix}\\ D_s = \begin{bmatrix} \vec x_1-\vec x_4 & \vec x_2-\vec x_4 & \vec x_3-\vec x_4 \end{bmatrix}\\ \begin{aligned} &\frac{\partial (D_s)_{kl}}{\partial \vec {\mathbf{x}_{ij}}} e_i\otimes e_j\otimes e_k\otimes e_l\\ &=\frac{\partial \vec {\mathbf{x}}_{kl}}{\partial \vec {\mathbf{x}}_{ij}} e_i\otimes e_j\otimes e_k\otimes e_l\\ &= \delta_{ik}\delta_{jl} e_i\otimes e_j\otimes e_k\otimes e_l \end{aligned} \] \(D_m^{-1}\)的分量表示为\(d_{mn}\)，\(P\)的分量表示为\(P_{rs}\)，则能量密度函数\(\Psi\)关于位置\(\vec {\mathbf{x}}\) 的梯度为：

\[ \begin{aligned} \frac{\partial \Psi}{\partial \vec {\mathbf{x}}} &= \frac{\partial F}{\partial \vec {\mathbf{x}}} : \frac{\partial \Psi}{\partial F}\\ &= : \frac{\partial F}{\partial \vec {\mathbf{x}}}: P \\ &= : (\frac{\partial D_s}{\partial \vec {\mathbf{x}}} D_m^{-1}): P \\ &= : (\delta_{ik}\delta_{jl} e_i\otimes e_j\otimes e_k\otimes e_l\cdot (d_{mn} e_m\otimes e_n)): P_{rs} e_r\otimes e_s \\ &= : (\delta_{ik}\delta_{jl}\delta_{lm} d_{mn} e_i\otimes e_j\otimes e_k\otimes e_n): P_{rs} e_r\otimes e_s \\ &= : (\delta_{ik} d_{jn} e_i\otimes e_j\otimes e_k\otimes e_n): P_{rs} e_r\otimes e_s \\ &= P_{kn} \delta_{ik} d_{jn} e_i\otimes e_j\\ &= P_{in} d_{jn} e_i\otimes e_j\\ &= P D_m^{-T} \end{aligned} \]

Matlab验证

下面的验证里将\(\vec {\mathbf{x}}\)展开为一个12维的向量。 \[ F=D_s D_m^{-1}\\ D_s = \begin{bmatrix} \vec x_1-\vec x_4 & \vec x_2-\vec x_4 & \vec x_3-\vec x_4 \end{bmatrix}\\ D_m = \begin{bmatrix} \vec x_{10}-\vec x_{40} & \vec x_{20}-\vec x_{40} & \vec x_{30}-\vec x_{40} \end{bmatrix}\\ \vec {\mathbf{x}} = \begin{bmatrix} \vec x_1 \\ \vec x_2 \\ \vec x_3 \\ \vec x_4 \end{bmatrix}\in \mathbb{R}^{12} \]

所以

\[ \frac{\partial F}{\partial \vec {\mathbf{x}}} = \frac{\partial D_s}{\partial \vec {\mathbf{x}}} D_m^{-1} \]

\(\frac{\partial D_s}{\partial \vec {\mathbf{x}}}\)是一个\(3\times 3\times 12\)的张量：

\[ \frac{\partial D_s}{\partial \vec {\mathbf{x}}} = \begin{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 1 & 0 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}\\ \begin{bmatrix} -1 & 0 & 0 \\ -1 & 0 & 0 \\ -1 & 0 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & -1 & 0 \\ 0 & -1 & 0 \\ 0 & -1 & 0 \end{bmatrix}\\ \begin{bmatrix} 0 & 0 & -1 \\ 0 & 0 & -1 \\ 0 & 0 & -1 \end{bmatrix}\\ \end{bmatrix} \]

因为\(D_m, P\)都是固定的，所以把他们用符号表示就可以验证了：

syms p11 p12 p13 p21 p22 p23 p31 p32 p33 real
syms d11 d12 d13 d21 d22 d23 d31 d32 d33 real

DmInv = [d11, d12, d13; d21, d22, d23; d31, d32, d33];
dDsdx = sym(zeros(3, 3, 12));

% 填充前九个矩阵，每个矩阵在对应位置置为1
for k = 1:9
    i = mod(k-1, 3) + 1;  % 计算行索引
    j = floor((k-1)/3) + 1;  % 计算列索引
    dDsdx(i, j, k) = 1;
end

% 填充最后三个矩阵
dDsdx(:,:,10) = [-1 -1 -1; 0 0 0; 0 0 0];
dDsdx(:,:,11) = [0 0 0; -1 -1 -1; 0 0 0];
dDsdx(:,:,12) = [0 0 0; 0 0 0; -1 -1 -1];

P = [p11, p12, p13; p21, p22, p23; p31, p32, p33];

dPdx = sym(zeros(12, 1));

% 计算double contraction
for k = 1:12
    % 算出vec(dFdx)前右乘DmInv
    dFdx = dDsdx(:,:,k) * DmInv;
    dFdx = dFdx(:);
    dPdx(k) = dFdx' * P(:);
end
dPdx = reshape(dPdx,3, 4);
dPdx(:,1:3) - P * DmInv'

最后发现dPdx的前三列和P * DmInv'是一样的，也就验证了这个公式。