# This is just a modified version of the MPM99 example script from taichi's simulation examples. import taichi as ti ti.init(arch=ti.cpu) dim = 2 quality = 1 # Use a larger value for higher-res simulations n_particles, n_grid = 9000 * 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 = 0.1e4, 0.2 # Young's modulus and Poisson's ratio # E, nu = 0.1e5, 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(dim, dtype=float, shape=n_particles) # position v = ti.Vector.field(dim, dtype=float, shape=n_particles) # velocity C = ti.Matrix.field(dim, dim, dtype=float, shape=n_particles) # affine velocity field F = ti.Matrix.field(dim, dim, 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) # plastic deformation grid_v = ti.Vector.field(dim, 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.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] # F[p]: deformation gradient update F[p] = (ti.Matrix.identity(float, 2) + dt * C[p]) @ F[p] # h: Hardening coefficient: snow gets harder when compressed h = ti.exp(10 * (1.0 - Jp[p])) if material[p] == 1: # jelly, make it softer h = 0.08 E, nu = 1000, 5 # Young's modulus and Poisson's ratio mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ( (1 + nu) * (1 - 2 * nu)) # Lame parameters if material[p] == 3: #stiff version h = 6 E, nu = 1e9, 0.3 mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ( (1 + nu) * (1 - 2 * nu)) # Lame parameters mu, la = mu_0 * h, lambda_0 * h if material[p] == 0: # liquid mu = 0 if material[p] == 4: #viscous mu = 0.1 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 = ti.min(ti.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 if material[p] == 0:#liquid # Reset deformation gradient to avoid numerical instability F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) # elif material[p] == 1: # F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) elif material[p] == 4: F[p] = U @ sig @ V.transpose() # F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J) elif material[p] == 2:#2 is snow # Reconstruct elastic deformation gradient after plasticity F[p] = U @ sig @ V.transpose() 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] # Loop over 3x3 grid node neighborhood for i, j in ti.static(ti.ndrange(3, 3)): 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 for i, j in grid_m: if grid_m[i, j] > 0: # No need for epsilon here grid_v[i, j] = \ (1 / grid_m[i, j]) * grid_v[i, j] # Momentum to velocity grid_v[i, j][1] -= dt * 50 # gravity 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 if j > n_grid - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0 if i == 0: grid_v[i,j][0] = 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] 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 group_size = n_particles // 5 @ti.kernel def initialize(): for i in range(n_particles): # x[i] = [ # ti.random() * 0.15 +0.1+ 0.15 * (i // group_size), # ti.random() * 0.15 + 0.05 + 0.32 * (1) # ] x[i] = [ ti.random() * 0.15 + 0.15 * (i // group_size), ti.random() * 0.15 + 0.05 + 0.32 * (1) ] if (i//group_size) == 4: material[i] = 4#i // group_size # 0: fluid 1: jelly 2: snow else: material[i] = 3 v[i] = ti.Matrix([0, 0]) F[i] = ti.Matrix([[1, 0], [0, 1]]) Jp[i] = 1 @ti.kernel def initgrid(): L = 0.5 W = 0.15 for i in range(n_particles): x[i] = [ti.random() * 0.5, ti.random()*0.15+0.5] #random version. material[i] = 3 v[i] = ti.Matrix([0, 0]) F[i] = ti.Matrix([[1, 0], [0, 1]]) Jp[i] = 1 def main(): initialize() # initgrid() gui = ti.GUI("Taichi MLS-MPM-99", res=512, background_color=0x112F41) while not gui.get_event(ti.GUI.ESCAPE, ti.GUI.EXIT): for s in range(int(2e-3 // dt)): substep() gui.circles(x.to_numpy(), radius=1.5, palette=[0x068587, 0xED553B, 0xEEEEF0,0xD2B4DE,0xEEEEF0], palette_indices=material) gui.get_event() if gui.is_pressed('r'): print('resetting') initialize() # initgrid() # Change to gui.show(f'{frame:06d}.png') to write images to disk gui.show() if __name__ == '__main__': main()