类似刚体球体在液体中流动，胡老师说使用相对速度边界条件，想仔细了解一下其原理。请问大家有demo或者相关学习资料吗？
Coulomb friction
需要用库伦摩擦模型。
我写了个刚体球和弹塑性体碰撞的程序，也许可以参考下，用鼠标控制球的运动：
（主要看 friction_coeff = 0.4 这行上下文的部分。）
import taichi as ti
import numpy as np
import random
ti.init(arch=ti.gpu) # Try to run on GPU
quality = 1 # Use a larger value for higher-res simulations
n_grid = (40 * quality, 80 * quality)
dx, inv_dx = 1 / n_grid[1], n_grid[1]
dt = 1e-4 / quality
p_vol, p_rho = (dx * 0.5)**2, 1
p_mass = p_vol * p_rho
E, nu = 0.15e4, 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
max_num_particles = 1024 * 16
dim = 2
x = ti.Vector.field(dim, dtype=float, shape=max_num_particles) # position
v = ti.Vector.field(dim, dtype=float, shape=max_num_particles) # velocity
C = ti.Matrix.field(dim, dim, dtype=float, shape=max_num_particles) # affine velocity field
F = ti.Matrix.field(dim, dim, dtype=float, shape=max_num_particles) # deformation gradient
material = ti.field(dtype=int, shape=max_num_particles) # material id
Jp = ti.field(dtype=float, shape=max_num_particles) # plastic deformation
# ti.root.dense(ti.i, max_num_particles).place(x, v, C, F, material, Jp)
cur_num_particles = ti.field(ti.i32, shape=())
collider_radius = 0.1
grid_v = ti.Vector.field(dim, dtype=float,
shape=n_grid) # grid node momentum/velocity
grid_m = ti.field(dtype=float, shape=n_grid) # grid node mass
@ti.kernel
def substep(collider_x: float, collider_y: float, collider_vx: float, collider_vy: float):
for i, j in grid_m:
grid_v[i, j] = [0, 0]
grid_m[i, j] = 0
# p2g
cur_num = cur_num_particles[None]
for p in range(cur_num):
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] = (ti.Matrix.identity(float, 2) +
dt * C[p]) @ F[p] # deformation gradient update
h = ti.exp(
10 *
(1.0 -
Jp[p])) # Hardening coefficient: snow gets harder when compressed
h = min(h, 15.0)
if material[p] >= 2: # jelly, make it softer
h = 0.3
mu, la = mu_0 * h, lambda_0 * h
if material[p] == 0: # liquid
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] == 1: # 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
if material[
p] == 0: # Reset deformation gradient to avoid numerical instability
F[p] = ti.Matrix.identity(float, 2) * ti.sqrt(J)
elif material[p] == 1:
F[p] = U @ sig @ V.transpose(
) # Reconstruct elastic deformation gradient after plasticity
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 node neighborhood
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
grid_pos = ti.Vector([i, j]) * dx
collider_center = ti.Vector([collider_x, collider_y])
collider_offset = grid_pos - collider_center
if collider_offset.norm() < collider_radius:
normal = collider_offset.normalized(1e-8)
v_co = ti.Vector([collider_vx, collider_vy]) # collider
v_rel = grid_v[i, j] - v_co
v_n = v_rel.dot(normal)
v_tan = v_rel - v_n * normal
v_tan_len = v_tan.norm()
friction_coeff = 0.4
if v_n < 0:
v_tan_normalized = v_tan / (v_tan_len + 1e-5)
v_rel = max(v_tan_len + friction_coeff * v_n, 0) * v_tan_normalized
grid_v[i, j] = v_rel + v_co
if i < 3 and grid_v[i, j][0] < 0:
grid_v[i, j][0] = 0 # Boundary conditions
if i > n_grid[0] - 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] = ti.Vector([0, 0])
if j > n_grid[1] - 3 and grid_v[i, j][1] > 0: grid_v[i, j][1] = 0
# g2p
for p in range(cur_num):
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
num_per_tetromino_square = 128
num_per_tetromino = num_per_tetromino_square * 4
staging_tetromino_x = ti.Vector.field(dim,
dtype=float,
shape=num_per_tetromino)
class StagingTetromino(object):
def __init__(self, x):
"""
Stores the info of the staging tetromino.
Args
x: a taichi field. Note this class just keeps |x| as a reference, and
is not a taichi data_oriented class.
"""
self._x_field = x
self.material_idx = 0
self._offsets = np.array([
[[0, -1], [1, 0], [0, -2]],
[[1, 1], [-1, 0], [1, 0]],
[[0, -1], [-1, 0], [0, -2]],
[[0, 1], [1, 0], [1, -1]],
[[1, 0], [2, 0], [-1, 0]],
[[0, 1], [1, 1], [1, 0]],
[[-1, 0], [1, 0], [0, 1]],
])
self._x_np_canonical = None
def regenerate(self, mat, kind):
self.material_idx = mat
shape = (num_per_tetromino, dim)
x = np.zeros(shape=shape, dtype=np.float32)
for i in range(1, 4):
# is there a more np-idiomatic way?
begin = i * num_per_tetromino_square
x[begin:(begin +
num_per_tetromino_square)] = self._offsets[kind, (i - 1)]
x += np.random.rand(*shape)
x *= 0.05
self._x_np_canonical = x
def update_center(self, center):
self._x_field.from_numpy(np.clip(self._x_np_canonical + center, 0, 1))
def rotate(self):
theta = np.radians(90)
c, s = np.cos(theta), np.sin(theta)
m = np.array([[c, -s], [s, c]], dtype=np.float32)
x = m @ self._x_np_canonical.T
self._x_np_canonical = x.T
staging_tetromino = StagingTetromino(staging_tetromino_x)
@ti.kernel
def drop_staging_tetromino(mat: int):
base = cur_num_particles[None]
for i in staging_tetromino_x:
bi = base + i
x[bi] = staging_tetromino_x[i]
material[bi] = mat
v[bi] = ti.Matrix([0, -2])
F[bi] = ti.Matrix([[1, 0], [0, 1]])
Jp[bi] = 1
cur_num_particles[None] += num_per_tetromino
def compute_tetromino_center(mouse):
# TODO: better clipping...
x = max(min(mouse[0], 0.45), 0.05)
return np.array([x, 0.8], dtype=np.float32)
def main():
gui = ti.GUI("Taichi Snow Effects",
res=(384, 768),
background_color=0x112F41)
collider_center = np.array([0.5, 0.5])
def gen_mat_and_kind():
material_id = 1
return material_id, random.randint(0, 6)
staging_tetromino.regenerate(*gen_mat_and_kind())
last_action_frame = -1e10
last_mouse = gui.get_cursor_pos()
for f in range(100000):
if gui.get_event(ti.GUI.PRESS):
ev_key = gui.event.key
if ev_key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: break
elif ev_key == ti.GUI.SPACE:
if cur_num_particles[
None] + num_per_tetromino < max_num_particles:
drop_staging_tetromino(staging_tetromino.material_idx)
print('# particles =', cur_num_particles[None])
staging_tetromino.regenerate(*gen_mat_and_kind())
last_action_frame = f
elif ev_key == 'r':
staging_tetromino.rotate()
mouse = gui.get_cursor_pos()
mouse = (mouse[0] * 0.5, mouse[1])
staging_tetromino.update_center(compute_tetromino_center(mouse))
collider_center = np.array(mouse)
if f == 0:
vx, vy = 0, 0
else:
vx = (mouse[0] - last_mouse[0]) / 2e-3
vy = (mouse[1] - last_mouse[1]) / 2e-3
last_mouse = mouse
for s in range(int(2e-3 // dt)):
mx, my = collider_center
substep(mx, my, vx, vy)
colors = np.array([
0xA6B5F7, 0xEEEEF0, 0xED553B, 0x3255A7, 0x6D35CB, 0xFE2E44,
0x26A5A7, 0xEDE53B
],
dtype=np.uint32)
particle_radius = 2.3
gui.circles(x.to_numpy() * [[2, 1]],
radius=particle_radius,
color=colors[material.to_numpy()])
if last_action_frame + 40 < f:
gui.circles(staging_tetromino_x.to_numpy() * [[2, 1]],
radius=particle_radius,
color=int(colors[staging_tetromino.material_idx]))
gui.circle(pos=collider_center * [2, 1], radius=collider_radius * 768)
gui.text('Taichi Snow Effects', (0.07, 0.97), font_size=20)
padding = 0.025
segments = 20
step = (1 - padding * 4) / (segments - 0.5) / 2
gui.line(begin=(padding * 2, padding),
end=(1 - padding * 2, padding),
radius=2)
gui.line(begin=(padding * 2, 1 - padding),
end=(1 - padding * 2, 1 - padding),
radius=2)
gui.line(begin=(padding * 2, padding),
end=(padding * 2, 1 - padding),
radius=2)
gui.line(begin=(1 - padding * 2, padding),
end=(1 - padding * 2, 1 - padding),
radius=2)
gui.show()
if __name__ == '__main__':
main()
2 Likes