-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmlsmpm_fluid1.py
116 lines (105 loc) · 4.59 KB
/
mlsmpm_fluid1.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
import taichi as ti
import numpy as np
ti.init(arch=ti.gpu)
n_particles, n_grid, n_lines = 9000, 128, 0
dx, inv_dx = 1 / n_grid, float(n_grid)
dt = 1e-4
p_vol, p_rho = (dx * 0.5)**2, 1
p_mass = p_vol * p_rho # particle mass
E, nu = 0.1e4, 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) # particle position
v = ti.Vector.field(2, dtype=float, shape=n_particles) # particle 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: 0,1,2: fluid
Jp = ti.field(dtype=float, shape=n_particles) # 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
@ti.func
def P2G(): # Particle to grid (P2G)
for p in x:
base = (x[p] * inv_dx - 0.5).cast(int)
fx = x[p] * inv_dx - base.cast(float)
# Bspline
w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2]
F[p] = (ti.Matrix.identity(float, 2) + dt * C[p]) @ F[p] # deformation gradient update
h = max(0.1, min(5, ti.exp(10 * (1.0 - Jp[p])))) # Hardening coefficient
mu, la = mu_0 * h, lambda_0 * h
mu = 0.0
U, sig, V = ti.svd(F[p])
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
sig[d, d] = new_sig
J *= new_sig
F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) # Reset deformation gradient
stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose() + ti.Matrix.identity(float, 2) * la * J * (J - 1)
stress = (-dt * p_vol * 4 * inv_dx * inv_dx) * stress
affine = stress + p_mass * C[p]
for i, j in ti.static(ti.ndrange(3, 3)): # Loop over 3x3 grid
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)
grid_m[base + offset] += weight * p_mass
@ti.func
def G2P(): # grid to particle (G2P)
for p in x:
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]
v_pic = ti.Vector.zero(float, 2)
C_pic = ti.Matrix.zero(float, 2, 2)
for i, j in ti.static(ti.ndrange(3, 3)): # loop over 3x3 grid
dpos = ti.Vector([i, j]).cast(float) - fx
g_v = grid_v[base + ti.Vector([i, j])]
weight = w[i][0] * w[j][1]
v_pic += weight * g_v
C_pic += 4 * inv_dx * weight * g_v.outer_product(dpos)
v[p], C[p] = v_pic, C_pic
x[p] += dt * v[p] # advection
@ti.func
def update_gridv():
for i, j in grid_m:
if grid_m[i, j] > 0:
grid_v[i, j] = (1 / grid_m[i, j]) * grid_v[i, j] # momentum to velocity
grid_v[i, j][1] -= dt * 250 # gravity
@ti.func
def enforce_boundary():
for i, j in grid_m:
# grid boundary
if i < 3 and grid_v[i, j][0] < 0: grid_v[i, j][0] = 0
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
if j > n_grid - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0
@ti.kernel
def substep():
for i, j in grid_m:
grid_v[i, j] = [0, 0]
grid_m[i, j] = 0
P2G()
update_gridv()
enforce_boundary()
G2P()
@ti.kernel
def initialize():
for i in range(n_particles):
x[i] = [ti.random() * 0.25 + 0.05, ti.random() * 0.25 + 0.05 ]
material[i] = 0 # 0,1,2: fluid
v[i] = ti.Matrix([0, 0])
F[i] = ti.Matrix([[1, 0], [0, 1]])
Jp[i] = 1
initialize()
gui = ti.GUI("Fluid MLS-MPM", res=512, background_color=0x000000)
while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT):
for s in range(int(2e-3 // dt)):
substep()
colors = np.array([0x36eeff, 0xfca13f, 0xEEEEF0], dtype=np.uint32)
gui.circles(x.to_numpy(), radius=1.5, color=colors[material.to_numpy()])
# gui.lines(lines1.to_numpy(), lines2.to_numpy(), color = 0x106b06, radius = 1.5)
gui.show() # only show gui
# gui.show(f'{gui.frame:06d}.png') # save image screenshots in current directory