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

轉載：  http://www.javashuo.com/article/p-fmwhagwr-md.html

彈簧質點模型的求解方法包括顯式歐拉積分和隱式歐拉積分等方法，其中顯式歐拉積分求解快速，但積分步長小，兩個可視幀之間須要屢次積分，而隱式歐拉積分則須要求解線性方程組，但其穩定性好，可以取較大的積分步長。[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

相關：

彈簧質點系統（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.