# 代码预览

## 线性系统求解部分

### 变量定义

``````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

``````
3 Likes

f_ij 的 变量 是 p_i , p_j (点i和点j的位置，二维向量)；
f_ij 的 值 是 力向量(二维向量)；

PS：这段代码写错了。