-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmpm128.py
161 lines (151 loc) · 7.48 KB
/
mpm128.py
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
import taichi as ti
ti.init(arch=ti.gpu) # Try to run on GPU
quality = 1 # Use a larger value for higher-res simulations
n_particles, n_grid = 300 * quality**2, 128 * quality
dx, inv_dx = 1 / n_grid, float(n_grid)
dt = 1e-4 / quality
p_vol, p_rho = (dx * 0.5) ** 2, 1
p_mass = p_vol * p_rho
E, nu = 5e3, 0.2 # Young's modulus and Poisson's ratio
mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu)) # Lame parameters 也就是超弹性模型的两个系数
x = ti.Vector.field(2, dtype=float, shape=n_particles) # position
v = ti.Vector.field(2, dtype=float, shape=n_particles) # velocity
C = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # affine velocity field
F = ti.Matrix.field(2, 2, dtype=float, shape=n_particles) # deformation gradient
material = ti.field(dtype=int, shape=n_particles) # material id
Jp = ti.field(dtype=float, shape=n_particles) # det of plastic deformation 通过直接维护应变行列式防止精度爆炸
grid_v = ti.Vector.field(2, dtype=float, shape=(n_grid, n_grid)) # grid node momentum/velocity
grid_m = ti.field(dtype=float, shape=(n_grid, n_grid)) # grid node mass
gravity = ti.Vector.field(2, dtype=float, shape=())
attractor_strength = ti.field(dtype=float, shape=())
attractor_pos = ti.Vector.field(2, dtype=float, shape=())
@ti.kernel
def substep():
for i, j in grid_m:
grid_v[i, j] = [0, 0]
grid_m[i, j] = 0
for p in x: # Particle state update and scatter to grid (P2G)
base = (x[p] * inv_dx - 0.5).cast(int)
fx = x[p] * inv_dx - base.cast(float)
# Quadratic kernels [http://mpm.graphics Eqn. 123, with x=fx, fx-1,fx-2]
w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2]
# deformation gradient update
F[p] = (ti.Matrix.identity(float, 2) + dt * C[p]) @ F[p] # 在这里C是速度对位移的梯度，在dt内积累的形变
# Hardening coefficient: snow gets harder when compressed
h = ti.max(0.1, ti.min(5, ti.exp(10 * (1.0 - Jp[p])))) # 这里雪最多被压缩到使 lame parameters 到 0.1 的境地
if material[p] == 1: # jelly, make it softer
h = 0.3
mu, la = mu_0 * h, lambda_0 * h # h 直接作用于 雪
if material[p] == 0: # liquid
mu = 0.0 # 在这里 liquid 根本不可压缩，因此直接将其mu设置为0
U, sig, V = ti.svd(F[p]) # 用svd 将旋转与拉伸进行分解
J = 1.0
for d in ti.static(range(2)):
new_sig = sig[d, d]
if material[p] == 2: # Snow
new_sig = min(max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3) # Plasticity
Jp[p] *= sig[d, d] / new_sig # 对det 进行更新
sig[d, d] = new_sig
J *= new_sig # 用trace 来更新J ，这里的J 是用来计算受力的，雪不需要这个
if material[p] == 0:
# Reset deformation gradient to avoid numerical instability
F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) # 对于水将旋转量消除因为这一部分主要计算压力，旋转量在P2G的affine matrix 中体现
elif material[p] == 2:
# Reconstruct elastic deformation gradient after plasticity
F[p] = U @ sig @ V.transpose() # 对于雪要把他的deformation进行reset 塑性形变要求我们去随时清除积累应变
stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose() + ti.Matrix.identity(float, 2) * la * J * (
J - 1
) # stress 是对势能的求导
stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress # 计算pressure of grid
affine = stress + p_mass * C[p] # 通过pressure 来更新affine
for i, j in ti.static(ti.ndrange(3, 3)):
# Loop over 3x3 grid node neighborhood 这一步是p2g
offset = ti.Vector([i, j])
dpos = (offset.cast(float) - fx) * dx
weight = w[i][0] * w[j][1]
grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos) # 角动量现已加入套餐，因为affine记录了各个点的位置
grid_m[base + offset] += weight * p_mass
for i, j in grid_m:
if grid_m[i, j] > 0: # No need for epsilon here
# Momentum to velocity
grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j]
grid_v[i, j] += dt * gravity[None] * 30 # gravity
dist = attractor_pos[None] - dx * ti.Vector([i, j])
grid_v[i, j] += dist / (0.01 + dist.norm()) * attractor_strength[None] * dt * 100
if i < 3 and grid_v[i, j][0] < 0:
grid_v[i, j][0] = 0 # Boundary conditions
if i > n_grid - 3 and grid_v[i, j][0] > 0:
grid_v[i, j][0] = 0
if j < 3 and grid_v[i, j][1] < 0:
grid_v[i, j][1] = 0 # 如果是底边就把y方向速度设置为0
if j > n_grid - 3 and grid_v[i, j][1] > 0:
grid_v[i, j][1] = 0 # 如果是顶部同理
for p in x: # grid to particle (G2P)
base = (x[p] * inv_dx - 0.5).cast(int)
fx = x[p] * inv_dx - base.cast(float)
w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1.0) ** 2, 0.5 * (fx - 0.5) ** 2] # B-spline插值函数
new_v = ti.Vector.zero(float, 2)
new_C = ti.Matrix.zero(float, 2, 2)
for i, j in ti.static(ti.ndrange(3, 3)):
# loop over 3x3 grid node neighborhood
dpos = ti.Vector([i, j]).cast(float) - fx
g_v = grid_v[base + ti.Vector([i, j])]
weight = w[i][0] * w[j][1]
new_v += weight * g_v
new_C += 4 * inv_dx * weight * g_v.outer_product(dpos)
v[p], C[p] = new_v, new_C
x[p] += dt * v[p] # advection
@ti.kernel
def reset():
group_size = n_particles // 3
for i in range(n_particles):
x[i] = [
ti.random() * 0.2 + 0.3 + 0.10 * (i // group_size),
ti.random() * 0.2 + 0.05 + 0.32 * (i // group_size),
]
material[i] = i // group_size # 0: fluid 1: jelly 2: snow
v[i] = [0, 0]
F[i] = ti.Matrix([[1, 0], [0, 1]])
Jp[i] = 1
C[i] = ti.Matrix.zero(float, 2, 2)
print("[Hint] Use WSAD/arrow keys to control gravity. Use left/right mouse buttons to attract/repel. Press R to reset.")
gui = ti.GUI("Taichi MLS-MPM-128", res=1024, background_color=0x112F41)
reset()
gravity[None] = [0, -2]
for frame in range(20000):
if gui.get_event(ti.GUI.PRESS):
if gui.event.key == "r":
reset()
elif gui.event.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]:
break
if gui.event is not None:
gravity[None] = [0, 0] # if had any event
if gui.is_pressed(ti.GUI.LEFT, "a"):
gravity[None][0] = -1
if gui.is_pressed(ti.GUI.RIGHT, "d"):
gravity[None][0] = 1
if gui.is_pressed(ti.GUI.UP, "w"):
gravity[None][1] = 1
if gui.is_pressed(ti.GUI.DOWN, "s"):
gravity[None][1] = -1
mouse = gui.get_cursor_pos()
gui.circle((mouse[0], mouse[1]), color=0x336699, radius=15)
attractor_pos[None] = [mouse[0], mouse[1]]
attractor_strength[None] = 0
if gui.is_pressed(ti.GUI.LMB):
attractor_strength[None] = 1
if gui.is_pressed(ti.GUI.RMB):
attractor_strength[None] = -1
for s in range(int(2e-3 // dt)):
substep()
gui.circles(
x.to_numpy(),
radius=1,
palette=[0x068587, 0xED553B, 0xEEEEF0],
palette_indices=material,
)
# if frame%100==0:
# print(F[8000])
# print(Jp[8000])
# Change to gui.show(f'{frame:06d}.png') to write images to disk
gui.show()