position based dynamics

应读者需求，今天终于抽出时间把这篇文章整理一下，不过只把其中的一些公式进行了推导，如果需要细看，仍需阅读论文

Position Based Dynamics文章来源

首先对比一下PDF与FEM之间的关系 

PBD是预估点的下一位置，然后通过constrain把它project到最终位置，并用这个位置去更新速度

FEM是算出每一帧中每个点的内力和外力，然后用合力去更新点的速度，进而再根据速度去更新位置，然后循环

3.3 Constraint Projection

4.4 Cloth Balloons

Appendix A

其中

最后上张撕裂效果图：

h=1/30;
%gravity=[0,-10,0];

X0=read_obj('cloth.obj'); % seahorse_extended.obj');
v0=X0.xyz';  %坐标
originalv0=v0;
pre_v0=v0;
tri=X0.tri'; %三角形顺序索引
velocity=v0;
gravity=zeros(X0.node_num,3);
gravity(:,2)=-10;
velocity(:,1)=0;velocity(:,2)=0;velocity(:,3)=0;  %速度初始化
%controlpoint = [1,3,7,9];
%controlpoint = [1,11,111,121];
controlpoint = [1,217];  %固定点的索引号.
%velocity(controlpoint(1,:),:)=-1000;

win_min = [999 999 999];
win_max = [-999 -999 -999];

timestep=200;  %total step
xAxis=[1:timestep];
for index=1:timestep
    velocity=velocity+gravity*h;
    v0=pre_v0+h*velocity;    %PBD中的位置预测
    for xx=1:3  %迭代次数
        v0=distanceconstraint(v0,X0,controlpoint);  %距离约束
    end
     %如果撞到地面的话（-10），按照上一时刻的位置，然后速度反向，弹回去,布料这个例子中并未用到与地面的碰撞%
    for i=1:X0.node_num
        if v0(i,2)<-20
            v0(i,2)=pre_v0(i,2)-velocity(i,2)*h;
        end
        %if v0(i,2)<-10
        %    v0(i,2)=-10;
       % end
    end
   
    %----------------------------------------------------------------%
    v0(controlpoint(1,:),:)= pre_v0(controlpoint(1,:),:);  %更新速度与位置
    velocity=(v0-pre_v0)/h;
    pre_v0=v0;
    center=mean(v0,1);

    trimesh(tri,v0(:,1),v0(:,2),v0(:,3),'EdgeColor','b'); 
    axis equal
    axis manual
    cur_V = v0(:,1:3);
    

    
    % scale window to fit
    win_min = min([win_min; cur_V]);
    win_max = max([win_max; cur_V]);
    axis(reshape([win_min;win_max],1,2*size(cur_V,2)))
    
    pause(0.01);
end
grid on;
%figure, trimesh(tri,v0(:,1),v0(:,2),v0(:,3),'EdgeColor','b'); axis equal
%hold on, plot3(v0(controlpoint(1,:),1),v0(controlpoint(1,:),2),v0(controlpoint(1,:),3),'*r');

function P=distanceconstraint(P,X0,controlpoint)
for i=1:X0.face_num
    idx1=X0.tri(1,i);  %三角形的三个索引
    idx2=X0.tri(2,i);
    idx3=X0.tri(3,i);
    d1=norm(X0.xyz(:,idx1)-X0.xyz(:,idx2));  %原始两点间的距离
    c_d1=norm(P(idx1,:)-P(idx2,:));%现在两点间的距离
    d2=norm(X0.xyz(:,idx2)-X0.xyz(:,idx3));
    c_d2=norm(P(idx2,:)-P(idx3,:));
    d3=norm(X0.xyz(:,idx1)-X0.xyz(:,idx3));
    c_d3=norm(P(idx1,:)-P(idx3,:));
    if c_d1~=d1  %考察三角形第一条边， 如果现在的距离与原始的距离不相等，那么更新位置，就像弹簧一样
        test=1;%用来测试这个点是不是control point， 如果是control point 那么不更新位置
        test2=1;%同上
        for count=1:length(controlpoint)
            if idx1==controlpoint(count)
                test=0;  %这个点不是control point
            end
            if idx2==controlpoint(count)
                test2=0;
            end
        end
        if test==0
            deltap1=[0,0,0];
            deltap2=(c_d1-d1)*(P(idx1,:)-P(idx2,:))/(norm(P(idx1,:)-P(idx2,:)));
        else
            deltap1=-0.5*(c_d1-d1)*(P(idx1,:)-P(idx2,:))/(norm(P(idx1,:)-P(idx2,:))); %文章中的stifness是乘在这里的，我没乘，我按质点质量相等来算的，所以是0.5
        end
        if test2==0
            deltap2=[0,0,0];
            deltap1=-(c_d1-d1)*(P(idx1,:)-P(idx2,:))/(norm(P(idx1,:)-P(idx2,:)));
        else
            deltap2=0.5*(c_d1-d1)*(P(idx1,:)-P(idx2,:))/(norm(P(idx1,:)-P(idx2,:)));
        end
        P(idx1,:)=P(idx1,:)+deltap1;
        P(idx2,:)=P(idx2,:)+deltap2;
    end
    if c_d2~=d2  %考察第二条边
        test=1;
        test2=1;
        for count=1:length(controlpoint)
            if idx2==controlpoint(count)
                test=0;
            end
            if idx3==controlpoint(count)
                test2=0;
            end
        end
        if test==0
            deltap23=[0,0,0];
            deltap32=(c_d2-d2)*(P(idx2,:)-P(idx3,:))/(norm(P(idx2,:)-P(idx3,:)));
        else
            deltap23=-0.5*(c_d2-d2)*(P(idx2,:)-P(idx3,:))/(norm(P(idx2,:)-P(idx3,:)));
        end
        if test2==0
            deltap32=[0,0,0];
            deltap23=-(c_d2-d2)*(P(idx2,:)-P(idx3,:))/(norm(P(idx2,:)-P(idx3,:)));
        else
            deltap32=0.5*(c_d2-d2)*(P(idx2,:)-P(idx3,:))/(norm(P(idx2,:)-P(idx3,:)));
        end
        P(idx2,:)=P(idx2,:)+deltap23;
        P(idx3,:)=P(idx3,:)+deltap32;
    end
    if c_d3~=d3 %考察第三条边
        test=1;
        test2=1;
        for count=1:length(controlpoint)
            if idx1==controlpoint(count)
                test=0;
            end
            if idx3==controlpoint(count)
                test2=0;
            end
        end
        if test==0
            deltap13=[0,0,0];
            deltap31=(c_d3-d3)*(P(idx1,:)-P(idx3,:))/(norm(P(idx1,:)-P(idx3,:)));
        else
            deltap13=-0.5*(c_d3-d3)*(P(idx1,:)-P(idx3,:))/(norm(P(idx1,:)-P(idx3,:)));
        end
        if test2==0
            deltap31=[0,0,0];
            deltap13=-(c_d3-d3)*(P(idx1,:)-P(idx3,:))/(norm(P(idx1,:)-P(idx3,:)));
        else
            deltap31=0.5*(c_d3-d3)*(P(idx1,:)-P(idx3,:))/(norm(P(idx1,:)-P(idx3,:)));
        end
        P(idx1,:)=P(idx1,:)+deltap13;
        P(idx3,:)=P(idx3,:)+deltap31;
    end
end
end

MATLAB源代码

如果有什么疑问，欢迎留言交流