Homework1 弹簧质点系统 雅克比迭代

结果

在阻尼较小的情况下,质点会加速移动。
未命名

代码预览

线性系统求解部分

变量定义

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

完整代码 https://github.com/kphmd/games201/blob/master/HW1/HW1.py

2 Likes