100行代码实现PBD

本文使用matlab实现了一个简单的PBD场景:双摆系统,代码(加上注释)不超过100行。

物理仿真我一般都是用C++写的,但是一个简单的场景,再加上OpenGL渲染,就要不少代码。但是物理上,有些我们关注的就是一些核心的求解步骤。今天我试着用最简单的方式来说明PBD的基本原理和步骤,我不想关注渲染,那matlab就是一个很好的工具了。我只用了不到100行的代码就实现了一个简单的PBD场景:双摆系统,如下图所示。

PBD Spring
具体代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
clear;
clc;
% rest configuration
X0 = [0,2,0]';
X1 = [0,1,0]';
X2 = [0,0,0]';
% set inverse mass of spring particles
invM0 = 0.0;
invM1 = 1.0;
invM2 = 1.0;
% set velocities of spring particles
v0 = [0,0,0]';
v1 = [0,0,0]';
v2 = [0,0,0]';
% spring reset length
d1 = norm(X1-X0);
d2 = norm(X2-X1);

% deformed configuration
x0 = [0.0,2.0,0.0]';
x1 = [1.5,2.0,0.0]';
x2 = [3.0,2.0,0.0]';
%time step
h = 0.01;
% external force
f = [0,0,0]';
% gravity
g = [0,-9.8,0]';

% draw animation of springs
figure; grid on; axis equal; view(0,90);
axis([-2,2, -2,2, -2,2]);
for i = 1:10000
%semi-Euler method
v0 = v0 + h * invM0 * f;
v1 = v1 + h * invM1 * f + h * g;
v2 = v2 + h * invM2 * f + h * g;

% predicted positions
px0 = x0 + h * v0;
px1 = x1 + h * v1;
px2 = x2 + h * v2;

% constraint solving
for index = 1:5
[px0,px1] = solveDistanceConstraint(px0,px1,invM0,invM1,d1);
[px1,px2] = solveDistanceConstraint(px1,px2,invM1,invM2,d2);
end

% update velocities
v0 = (px0 - x0)/h;
v1 = (px1 - x1)/h;
v2 = (px2 - x2)/h;
% update positions
x0 = px0;
x1 = px1;
x2 = px2;
% draw springs
clf; hold on; grid on; axis equal; view(0, 90);
axis([-2,2, -1,3, -2,2]);
plotSpring(x0,x1,x2,'-ob');
drawnow();
end

function [px0,px1] = solveDistanceConstraint(px0,px1,invM0, invM1,d)
% solveDistanceConstraint

% constraint
c = norm(px0-px1) - d;

% compute gradient
g0 = (px1-px0)/norm(px1-px0);
g1 = -g0;
gradient = [g0;g1];

% inverse Mass matrix
v = [invM0,invM0,invM0,invM1,invM1,invM1];
invM = diag(v);

% common part: lagrangian multiplier
lag = (gradient' * invM * gradient) \ c;

deltaX0 = invM0 * lag * g0;
deltaX1 = invM1 * lag * g1;

px0 = px0 + deltaX0;
px1 = px1 + deltaX1;
end

function plotSpring(x0,x1,x2,color)
% plotSpring
M = zeros(3,3);
M(:,1) = x0;
M(:,2) = x1;
M(:,3) = x2;
order = [1,2,2,3];
plot3(M(1,order),M(2,order),M(3,order),color);
end

参考文献:

Position based dynamics

文章目录
  1. 1. 参考文献:
,