position based dynamics

position based dynamics implemented in matlab

function P=distanceconstraint（PX0controlpoint）
for i=1:X0.face_num
    idx1=X0.tri（1i）;  %三角形的三个索引
    idx2=X0.tri（2i）;
    idx3=X0.tri（3i）;
    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=[000];
            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=[000];
            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=[000];
            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=[000];
            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=[000];
            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=[000];
            deltap13=-（c_d3-d3）*（P（idx1:）-P（idx3:））/（norm（P（idx1:）-P（idx3:）））;
        else
            de

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件      33041  2015-01-27 21:49  PBD\cloth.obj

     文件       3325  2015-01-27 21:41  PBD\distanceconstraint.m

     文件       1711  2016-07-05 00:18  PBD\PBD.asv

     文件       1762  2016-07-05 00:33  PBD\PBD.m

     文件      30306  2013-03-29 16:34  PBD\read_obj.m

     文件       4378  2013-03-29 18:07  PBD\write_obj.m

     目录          0  2016-07-05 00:18  PBD

----------- ---------  ---------- -----  ----

                74523                    7