[Taichi入门]基于Taichi的简单经验刚体模拟

从零开始学Taichi的第一步,实现了一个基于经验的简单刚体系统。

由于仅包含球体,这个刚体系统实际上只有四条基本的经验法则,即:

  • 物体会受到重力的影响往下落
  • 物体和其它球体相撞会受到与撞击方向相反的作用力
  • 物体碰到墙壁会向相反方向反弹
  • 碰撞会使得力会按一定比例衰减

将这些经验法则抽象为数学公式,我们可以得到这样的结论:

考虑法则1,对于每一个子步,物体受到重力影响可以表示为:

考虑法则2,由于我们的模型里只有球体,因此当两物体相撞时,我们可以将两个球体的球心连线,这条连线就是球体的受力方向。同时,若两个球体距离过近,就说明其位置信息需要立刻纠正,因此受到的碰撞力度将更大。抽象为公式,则有:

k为控制碰撞强度的比例系数。

考虑法则3,当物体与墙壁碰撞时,直接按照速度与墙面的夹角变换即可。一种最简单的模拟方法是,考虑墙面始终与坐标轴垂直,当球体与某一方向的墙面碰撞时,直接使对应方向的速度矢量变为相反数即可。以球体和地面碰撞为例:

其中k为碰撞损耗能量的系数。

考虑法则4,使得物体的受力一直衰减即可。

结合上述公式,可以写出如下子步函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
@ti.kernel
def substep():
for i in range(n):
V[i] = V[i] + gravity * dt
if (abs(P[i].y + V[i].y * dt)) >= 1:
V[i].y = -V[i].y * 0.85
if (abs(P[i].x + V[i].x * dt)) >= Bound:
V[i].x = -V[i].x * 0.85
if (abs(P[i].z + V[i].z * dt)) >= Bound:
V[i].z = -V[i].z * 0.85
for j in range(n):
if (P[i] + V[i] * dt - P[j]).norm() <= r * 2:
V[i] = V[i] + (P[i] + V[i] * dt - P[j]) * 0.02
V[i] = V[i] * 0.99997
P[i] += V[i] * dt

完整代码如下:

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
import datetime
import time

import taichi as ti
import taichi.math as tm

# common scope
ti.init(arch=ti.gpu)
n = 1024
r = 0.01
dt = 1e-1 / n
substeps = int(1 / 60 // dt)

# physical scope
P = ti.Vector.field(3, dtype=float, shape=(n,))
V = ti.Vector.field(3, dtype=float, shape=(n,))
gravity = ti.Vector([0, -9.8, 0])
Bound = 1


# init
@ti.kernel
def initialize():
for i in range(n):
P[i] = [ti.random(float), ti.random(float), ti.random(float)]
V[i] = [ti.random(float) * 0.2, 0, ti.random(float) * 0.2]


# substep
@ti.kernel
def substep():
for i in range(n):
V[i] = V[i] + gravity * dt
if (abs(P[i].y + V[i].y * dt)) >= 1:
V[i].y = -V[i].y * 0.85
if (abs(P[i].x + V[i].x * dt)) >= Bound:
V[i].x = -V[i].x * 0.85
if (abs(P[i].z + V[i].z * dt)) >= Bound:
V[i].z = -V[i].z * 0.85
for j in range(n):
if (P[i] + V[i] * dt - P[j]).norm() <= r * 2:
V[i] = V[i] + (P[i] + V[i] * dt - P[j]) * 0.02
V[i] = V[i] * 0.99997
P[i] += V[i] * dt


# rendering
window = ti.ui.Window("Taichi Cloth Simulation on GGUI", (800, 800), show_window=True)

canvas = window.get_canvas()
canvas.set_background_color((0.3, 0.3, 0.3))
scene = ti.ui.Scene()
camera = ti.ui.Camera()

current_t = 0.0
initialize()
flag = 0

while window.running:
p1 = time.time()
if current_t > 8:
# Reset
initialize()
current_t = 0
flag = 0
# break

for i in range(substeps):
substep()
current_t += dt

camera.position(2, 1, 3)
camera.lookat(0.0, -0.5, 0)
scene.set_camera(camera)

scene.point_light(pos=(0, 1, 2), color=(1, 1, 1))
scene.ambient_light((0.4, 0.4, 0.4))

# Draw a smaller ball to avoid visual penetration
scene.particles(P, radius=r, color=(0.5, 0.42, 0.8))
canvas.scene(scene)
# filename = "output" + str(flag) + ".jpg"
# window.save_image(filename)
window.show()
flag += 1

p2 = time.time()
delta_time = (datetime.datetime.fromtimestamp(p2) - datetime.datetime.fromtimestamp(p1)).microseconds / 1000000.0

# print(delta_time, ' and ', 1. / 60)
if delta_time < 1. / 60:
time.sleep(1. / 60 - delta_time)

模拟效果如下: