结果
在阻尼较小的情况下，质点会加速移动。
代码预览
线性系统求解部分
变量定义
A = ti.Matrix(2, 2, dt=ti.f32, shape=(max_num_particles, max_num_particles))
b = ti.Vector(2, dt=ti.f32, shape=max_num_particles)
r = ti.Vector(2, dt=ti.f32, shape=())
dv = ti.Vector(2,dt=ti.f32, shape=max_num_particles)
new_dv = ti.Vector(2,dt=ti.f32, shape=max_num_particles)
雅克比迭代
@ti.func
def iterate(): # 更新dv，其中dv = v[i+1] - v[i]
n = num_particles[None]
for i in range(n):
r[None].x = b[i].x
r[None].y = b[i].y
for j in range(n):
if i != j:
r[None] -= A[i, j] @ dv[j]
new_dv[i].x = (r[None].x - A[i,i][0,1] * dv[i].y) / A[i, i][0,0]
new_dv[i].y = (r[None].y - A[i,i][1,0] * dv[i].x) / A[i, i][1,1]
for i in range(n):
dv[i] = new_dv[i]
@ti.func
def resi() -> ti.f32: # 计算残差
n = num_particles[None]
res = 0.0
for i in range(n):
r = b[i] * 1.0
for j in range(n):
r -= A[i, j] @ dv[j]
res += r.x ** 2 + r.y ** 2
return res
@ti.func
def solve_equation() :
for i in range(1000000) :
iterate()
if resi() < 1e-9 :
break
线性系统构造部分
力向量F & 力向量F对位置向量x的导数
@ti.func
def f_ij(i,j) : # 质点j对质点i的力向量F
x_ij = x[j] - x[i] # from i to j
x_d = x_ij.normalized()
x_n = x_ij.norm()
return spring_stiffness[None] * (x_n - rest_length[i, j]) * x_d
@ti.func
def df_ij(i,j) : # 质点j对质点i的力向量F的导数
x_ij = x[j] - x[i] # from i to j
x_d = x_ij.normalized()
x_n = x_ij.norm()
x_o = x_d.outer_product(x_d)
x_oi = x_o - [[1,0],[0,1]]
return - spring_stiffness[None] * ((x_n - rest_length[i,j]) / x_n * x_oi - x_o)
构造线性系统系数矩阵
@ti.kernel
def substep_jacobi():
# fill A b
n = num_particles[None]
dt2 = dt * dt
for i in range(n) :
for j in range(n) :
if i==j :
A[i,j] = [[1,0],[0,1]]
elif rest_length[i, j] != 0:
A[i,j] = - dt2 * df_ij(i,j)
else :
A[i,j] = [[0,0],[0,0]]
for i in range(n) :
b[i] = ti.Vector(gravity) * dt
for j in range(n) :
if i!=j and rest_length[i,j] != 0 :
b[i] += dt * f_ij(i,j) + dt2 * df_ij(i,j) @ v[j]
# solve delta_v
solve_equation()
# update v
for i in range(n):
v[i] *= ti.exp(-dt * damping[None]) # damping
v[i] += dv[i]
# Collide with ground
for i in range(n):
if x[i].y < bottom_y:
x[i].y = bottom_y
v[i].y = 0
if x[i].y > top_y:
x[i].y = top_y
v[i].y = 0
# Collide with wall
for i in range(n):
if x[i].x < left_x:
x[i].x = left_x
v[i].x = 0
if x[i].x > right_x:
x[i].x = right_x
v[i].x = 0
# Compute new position
for i in range(num_particles[None]):
x[i] += v[i] * dt