網格彈簧質點系統模擬（Spring-Mass System by Fast Method）附源碼

　　彈簧質點模型的求解方法包括顯式歐拉積分和隱式歐拉積分等方法，其中顯式歐拉積分求解快速，但積分步長小，兩個可視幀之間需要多次積分，而隱式歐拉積分則需要求解線性方程組，但其穩定性好，能夠取較大的積分步長。[Liu et al. 2007]文章提出了一種彈簧質點模型的求解方法，它將隱式歐拉積分方法轉變為求解最優化問題，並采用迭代分步優化的方法來達到最優解。相比隱式歐拉積分，該方法計算快速，並且精度在可接受范圍內。

　　彈簧質點模型的隱式表達方式如下：

（1）

（2）

其中：q_n和v_n分別代表t_n時刻質點的位置和速度，f(q_n)為t_n時刻質點所受到的力，M為質點的質量，h為步長。

　　利用式（1）我們可以得到：

（3）

（4）

　　將式（3）減式（4）並與式（2）結合得到：

（5）

　　記x = q_n+1，y = 2q_n – q_n-1，式（5）可以變化為：

（6）

　　式（6）的解其實對應於如下函數的臨界點：

（7）

　　於是彈簧質點模型問題可以變化為最優化問題min_x g(x)，即最小化函數g(x)。

　　函數E(x)中最重要的部分是彈簧勢能，根據Hooke定律，可以推導得到兩個質點間彈簧的勢能為：

（8）

其中：k為彈簧的彈性系數，r為彈簧的自然長度。

　　因此彈簧質點模型中彈簧的整體勢能也可以變化為最優化問題，即最小化如下函數：

（9）

其中：L = A·K·A^T，J = A·K，式中A∈R^m×s（m為質點數量，s為彈簧數量），並且A_i1,i=1，A_i2,i= -1，K∈R^m×m為對角矩陣，K_i,i = k_i。

　　如果考慮其他外力（如重力等），那么函數E(x)的表達式為：

（10）

其中：是所有彈簧為自然長度時的方向。

　　將函數E(x)的表達式（10）代入式（7），整理后得到最終的優化表達式：

（11）

　　對於上述優化問題，可以分兩步進行，將前一時刻的質點位置作為初始值x，首先固定x優化d，然后固定d優化x，然后重復上述迭代步驟直到滿足設定的迭代步數。

 

function [X, V] = spring_mass_fast(X0, V0, E, b, bc, R, h)
    % This code implements algorithm of the following paper:
    % "Fast Simulation of Mass-Spring Systems"

    m = size(X0,1); % vertex number
    s = size(E,1); % spring number
    
    if ~exist('R', 'var')
        R = normrow(X0(E(:,1),:) - X0(E(:,2),:));
    end
    
    damping = 0.02;
    drag = 1 - damping;
    stiffness = 1e1;
    K = stiffness*ones(s,1);
    mass = 0.01;
    M = diag(mass*ones(m,1));
    g = [0 0 -9.8];
    fext = repmat(mass*g, [m,1]);

    A = sparse(E,[1:s;1:s]',repmat([1,-1],s,1),m,s);
    
    L = A*diag(K)*A';
    J = A*diag(K);

    X = X0;
    iter = 0;
    max_iter = 10;
    while true
        % step1: Fix X and find D
        D = X(E(:,1),:) - X(E(:,2),:);
        D = bsxfun(@times, D, R./normrow(D));
        
        % step2: Fix D and find X
        X = solve_equation(M + h^2*L, h^2*(fext + J*D) + M*(X0 + V0*h), b, bc);
        
        iter = iter + 1;
        if iter == max_iter
            break;
        end
    end
    V = drag*(X - X0)/h;
end

 

本文為原創，轉載請注明出處：http://www.cnblogs.com/shushen。

 

 

相關： 

彈簧質點系統（Euler Integration）：http://www.cnblogs.com/shushen/p/5473264.html

彈簧質點系統（Verlet Integration）：http://www.cnblogs.com/shushen/p/5394431.html

參考文獻：

[1] Tiantian Liu, Adam W. Bargteil, James F. O'Brien, and Ladislav Kavan. 2013. Fast simulation of mass-spring systems. ACM Trans. Graph. 32, 6, Article 214 (November 2013), 7 pages.