Skip to content

Instantly share code, notes, and snippets.

@rexwangcc
Last active September 3, 2021 23:12
Show Gist options
  • Save rexwangcc/95be09722474aa3c66d2b59faea21266 to your computer and use it in GitHub Desktop.
Save rexwangcc/95be09722474aa3c66d2b59faea21266 to your computer and use it in GitHub Desktop.
A Moving Least Squares Material Point Method implementation with Taichi language. Originally created by @yuanming-hu with 88 lines of code.
import taichi as ti
ti.init(arch=ti.gpu)
n_particles = 8192
n_grid = 128
dx = 1 / n_grid
dt = 2e-4
p_rho = 1
p_vol = (dx * 0.5)**2
p_mass = p_vol * p_rho
gravity = 9.8
bound = 3
E = 400
x = ti.Vector.field(2, float, n_particles)
v = ti.Vector.field(2, float, n_particles)
C = ti.Matrix.field(2, 2, float, n_particles)
J = ti.field(float, n_particles)
grid_v = ti.Vector.field(2, float, (n_grid, n_grid))
grid_m = ti.field(float, (n_grid, n_grid))
@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:
Xp = x[p] / dx
base = int(Xp - 0.5)
fx = Xp - base
w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
stress = -dt * 4 * E * p_vol * (J[p] - 1) / dx**2
affine = ti.Matrix([[stress, 0], [0, stress]]) + p_mass * C[p]
for i, j in ti.static(ti.ndrange(3, 3)):
offset = ti.Vector([i, j])
dpos = (offset - fx) * dx
weight = w[i].x * w[j].y
grid_v[base + offset] += weight * (p_mass * v[p] + affine @ dpos)
grid_m[base + offset] += weight * p_mass
for i, j in grid_m:
if grid_m[i, j] > 0:
grid_v[i, j] /= grid_m[i, j]
grid_v[i, j].y -= dt * gravity
if i < bound and grid_v[i, j].x < 0:
grid_v[i, j].x = 0
if i > n_grid - bound and grid_v[i, j].x > 0:
grid_v[i, j].x = 0
if j < bound and grid_v[i, j].y < 0:
grid_v[i, j].y = 0
if j > n_grid - bound and grid_v[i, j].y > 0:
grid_v[i, j].y = 0
for p in x:
Xp = x[p] / dx
base = int(Xp - 0.5)
fx = Xp - base
w = [0.5 * (1.5 - fx)**2, 0.75 - (fx - 1)**2, 0.5 * (fx - 0.5)**2]
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)):
offset = ti.Vector([i, j])
dpos = (offset - fx) * dx
weight = w[i].x * w[j].y
g_v = grid_v[base + offset]
new_v += weight * g_v
new_C += 4 * weight * g_v.outer_product(dpos) / dx**2
v[p] = new_v
x[p] += dt * v[p]
J[p] *= 1 + dt * new_C.trace()
C[p] = new_C
@ti.kernel
def init():
for i in range(n_particles):
x[i] = [ti.random() * 0.4 + 0.2, ti.random() * 0.4 + 0.2]
v[i] = [0, -1]
J[i] = 1
init()
gui = ti.GUI('MPM88')
while gui.running and not gui.get_event(gui.ESCAPE):
for s in range(50):
substep()
gui.clear(0x112F41)
gui.circles(x.to_numpy(), radius=1, color=0x068587)
gui.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment