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

　　模擬物體變形最簡單的方法就是采用彈簧質點系統（Spring-Mass System），由於模型簡單並且實用，它已被廣泛應用於服飾、毛發以及彈性固體的動態模擬。對於三角網格而言，彈簧質點系統將網格中的頂點看作系統中的質點，而網格的邊則是連接這些質點的彈簧。這樣，彈簧質點系統模型就將物體簡化成由彈簧和質點組成的系統，並利用彈簧質點的運動規律來描述物體的彈性變形過程。

　　Verlet積分是求解牛頓運動方程的數值方法，原理簡單描述如下：首先將系統t+dt時刻的位置x(t+dt)以及系統t-dt時刻的位置x(t-dt)用泰勒公式展開：

 

　　上面兩式相加后得到：

 

　　進一步變化得到：

 

　　因此通過上式可以根據系統前兩時刻的狀態求解系統的當前狀態，這與”基於網格的波動方程模擬“一文中的求解過程有些類似。

　　為了真實模擬物體變形效果，需要對彈簧質點系統進行受力分析：1. 每個質點有自身重力的影響；2. 每個質點受到與它相連的彈簧彈力影響，彈簧彈力遵守胡克定律；3. 質點運動時受到與其速度成正比的阻尼約束；4. 質點會受到其他外力的影響，由於施加的外力在每個三角面片上有一個法向分量，我們只需對每個質點周圍三角片上的這些分量相加即可。

 

% constrains option
wind = false;
ball = true;
pins = false;

figure('Position', [400, 400, 400, 320]);
fh = drawMesh(V,F,'facecolor','y','edgecolor','none');
if ball
    center = [50 60 -80];
    radius = 40;
    drawSphere([center radius], 'facecolor','r', 'nPhi',96, 'nTheta',48);
end
if pins
    plot3([-10;110], [0;0], [0;0], 'k-', 'linewidth',2);
end
view([-30 20])
axis equal
axis off
axis([-10 100 -10 100 -110 0]);
camlight
lighting gouraud

set(gca, 'position', [0 0 1 1]);

% initial condition
x_pre = V;
x_cur = V;

% rest length
E = edges(F);
l0 = vectorNorm3d(V(E(:,1),:) - V(E(:,2),:)); 

nV = size(V,1);
draw_t = 0;
tic;
while true
    % spring force
    Fs = stiffness * (vectorNorm3d(x_cur(E(:,1),:) - x_cur(E(:,2),:)) - l0);
    dir = normalizeVector3d(x_cur(E(:,2),:) - x_cur(E(:,1),:));

    M1 = sparse(E, E, [Fs.*dir(:,1);-Fs.*dir(:,1)]);
    M2 = sparse(E, E, [Fs.*dir(:,2);-Fs.*dir(:,2)]);
    M3 = sparse(E, E, [Fs.*dir(:,3);-Fs.*dir(:,3)]);
    as = [diag(M1), diag(M2), diag(M3)] ./ m;

    % wind force
    aw = zeros(nV,3);
    if wind
        N = normalizeVector3d(normals(x_cur,F));
        Fw = N * wind_force(i/10)' .* wind_strength;

        M1 = sparse(F, F, repmat(Fw.*N(:,1),1,3));
        M2 = sparse(F, F, repmat(Fw.*N(:,2),1,3));
        M3 = sparse(F, F, repmat(Fw.*N(:,3),1,3));
        aw = [diag(M1), diag(M2), diag(M3)] ./ m;
    end

    % verlet integration with a simple damping model
    x_new = drag*(x_cur - x_pre) + x_cur + bsxfun(@plus, as+aw, g)*dt*dt;

    x_pre = x_cur;
    x_cur = x_new;

    % ball constrains
    if ball
        diff = bsxfun(@minus, x_cur, center);
        index = vectorNorm3d(diff) < radius+1;
        x_cur(index,:) = bsxfun(@plus, center, bsxfun(@times, normalizeVector3d(diff(index,:)), radius+1));
    end

    % pin constrains
    if pins
        x_pre(pin_idx,:) = V(pin_idx,:);
        x_cur(pin_idx,:) = V(pin_idx,:);
    end

    % updata figure
    if toc > 0.033
        set(fh, 'Vertices', x_cur);
        drawnow;
        tic;
    end
end

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