Force computation in FEM

证明1

首先，已知

\[ \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} \] 这里把\(\vec {\mathbf{x}}\)后三列写成一个3x3的矩阵。\(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} \]

证明2

\[ \frac{\partial F}{\partial x} = \frac{\partial D_s}{\partial x}D_m^{-1}\\ \text{vec}(\frac{\partial F}{\partial x}) = (\begin{bmatrix} -1 & -1 & -1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} Dm^{-1})^T \otimes I_{3\times 3}\\ \]

\[ \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 \\ &= \text{vec}(\frac{\partial F}{\partial \vec {\mathbf{x}}})^T\text{vec}(P) \\ &= (\begin{bmatrix} -1 & -1 & -1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} Dm^{-1}) \otimes I_{3\times 3}\cdot \text{vec}(P)\\ &= \text{vec}(P (\begin{bmatrix} -1 & -1 & -1 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix} Dm^{-1})^{T})\\ &= \text{vec}(P D_m^{-T} \begin{bmatrix} -1 & 1 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ -1 & 0 & 0 & 1 \end{bmatrix})\\ \end{aligned} \]

倒数第二个等式是因为Kronecker积的性质:

给定矩阵\(A,B,X\)：

\[ B^T \otimes A \cdot \text{vec}(X) = \text{vec}(AXB) \]

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}}}\)是一个\(12\times 3\times 3\)的张量。

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

syms p11 p12 p13 p21 p22 p23 p31 p32 p33 real
x = sym('x', [12, 1]);
P = sym('p', [3, 3]);

Ds = [x(4)-x(1) x(7)-x(1) x(10)-x(1); x(5)-x(2) x(8)-x(2) x(11)-x(2); x(6)-x(3) x(9)-x(3) x(12)-x(3)];
DmInv = sym('DmInv', [3, 3]);
PDsPx = sym(zeros(12, 3, 3));
for i = 1:3
    for j = 1:3
        for k = 1:12
            PDsPx(k, i, j) = diff(Ds(i, j), x(k));
        end
    end
end

dPdx = sym(zeros(12, 1));
dFdx = tensorDot(PDsPx, DmInv);

dPsidx = P * transpose([-1 -1 -1; 1 0 0; 0 1 0; 0 0 1]*DmInv);
dPsidx2 = reshape(doubleDotProd(dFdx, P),3,4);
disp(dPsidx2-dPsidx);

附录：张量计算的matlab代码

代码由claude.ai生成。

双点积

function C = doubleDotProd(A, B, contractionDims)
    % Double contraction of tensors A and B
    % Input:
    %   A: First tensor (n-dimensional array)
    %   B: Second tensor (m-dimensional array)
    %   contractionDims: [optional] 2x2 array specifying which dimensions to contract
    %                    Default is last two dimensions of A with first two of B
    % Output:
    %   C: Resulting tensor after double contraction
    
    % Get sizes of input tensors
    sizeA = size(A);
    sizeB = size(B);
    orderA = ndims(A);
    orderB = ndims(B);
    
    % If contractionDims not specified, use default
    if nargin < 3
        contractionDims = [orderA-1, orderA; 1, 2];
    end
    
    % Verify dimensions match for contraction
    dimA1 = sizeA(contractionDims(1,1));
    dimA2 = sizeA(contractionDims(1,2));
    dimB1 = sizeB(contractionDims(2,1));
    dimB2 = sizeB(contractionDims(2,2));
    
    if dimA1 ~= dimB1 || dimA2 ~= dimB2
        error('Dimensions must match for contraction');
    end
    
    % Special case: when both inputs are 2nd order tensors
    if orderA == 2 && orderB == 2
        C = sum(sum(A .* B));
        return;
    end
    
    % Reshape A and B for matrix multiplication
    remainingDimsA = setdiff(1:orderA, contractionDims(1,:));
    remainingDimsB = setdiff(1:orderB, contractionDims(2,:));
    
    % Calculate output dimensions
    newSizeA = [prod(sizeA(remainingDimsA)), dimA1*dimA2];
    newSizeB = [dimB1*dimB2, prod(sizeB(remainingDimsB))];
    
    % Permute and reshape tensors
    permA = [remainingDimsA, contractionDims(1,:)];
    permB = [contractionDims(2,:), remainingDimsB];
    
    A_reshaped = reshape(permute(A, permA), newSizeA);
    B_reshaped = reshape(permute(B, permB), newSizeB);
    
    % Perform contraction via matrix multiplication
    C_temp = A_reshaped * B_reshaped;
    
    % Handle different cases for output reshaping
    if isempty(remainingDimsA) && isempty(remainingDimsB)
        % Case 1: No remaining dimensions (scalar output)
        C = C_temp;
    elseif isempty(remainingDimsA)
        % Case 2: Only remaining dimensions from B
        if isscalar(remainingDimsB)
            C = C_temp;  % Return as vector directly
        else
            C = reshape(C_temp, sizeB(remainingDimsB));
        end
    elseif isempty(remainingDimsB)
        % Case 3: Only remaining dimensions from A
        if isscalar(remainingDimsA)
            C = C_temp;  % Return as vector directly
        else
            C = reshape(C_temp, sizeA(remainingDimsA));
        end
    else
        % Case 4: Remaining dimensions from both A and B
        finalSize = [sizeA(remainingDimsA), sizeB(remainingDimsB)];
        C = reshape(C_temp, finalSize);
    end
end

张量点积

function C = tensorDot(A, B, dims)
    % Compute the dot product between two tensors along specified dimensions
    % Input:
    %   A: First tensor (n-dimensional array)
    %   B: Second tensor (m-dimensional array)
    %   dims: [optional] 1x2 array specifying which dimensions to contract
    %         Default is last dimension of A with first dimension of B
    % Output:
    %   C: Resulting tensor after dot product
    
    % Get sizes of input tensors
    sizeA = size(A);
    sizeB = size(B);
    
    % If dims not specified, use default
    if nargin < 3
        dims = [ndims(A), 1];
    end
    
    % Verify dimensions match for contraction
    if sizeA(dims(1)) ~= sizeB(dims(2))
        error('Contracted dimensions must match: %d != %d', ...
              sizeA(dims(1)), sizeB(dims(2)));
    end
    
    % Special case: if both inputs are vectors
    if isVector(A) && isVector(B)
        C = A(:)' * B(:);
        return;
    end
    
    % Get remaining dimensions
    remainingDimsA = setdiff(1:ndims(A), dims(1));
    remainingDimsB = setdiff(1:ndims(B), dims(2));
    
    % Permute arrays to move contracted dimensions to the end for A
    % and to the beginning for B
    permA = [remainingDimsA, dims(1)];
    permB = [dims(2), remainingDimsB];
    
    A_perm = permute(A, permA);
    B_perm = permute(B, permB);
    
    % Reshape tensors into matrices
    sizeA_new = [prod(sizeA(remainingDimsA)), sizeA(dims(1))];
    sizeB_new = [sizeB(dims(2)), prod(sizeB(remainingDimsB))];
    
    A_reshaped = reshape(A_perm, sizeA_new);
    B_reshaped = reshape(B_perm, sizeB_new);
    
    % Perform matrix multiplication
    C_temp = A_reshaped * B_reshaped;
    
    % Reshape result back to appropriate dimensions
    finalSize = [sizeA(remainingDimsA), sizeB(remainingDimsB)];
    
    % If result is a scalar or vector, don't reshape
    if numel(finalSize) <= 1
        C = C_temp;
    else
        C = reshape(C_temp, finalSize);
    end
end

function result = isVector(A)
    % Helper function to check if input is a vector
    sz = size(A);
    result = sum(sz > 1) <= 1;
end