【物理模拟】超简单的shape matching模拟刚体运动

Github repo

https://github.com/chunleili/tiRigidBody

Demo

How to run

First, intall latest taichi by

pip install taichi

Then, run with

python tiRigidBody

Rember to press SPACE to start the simulation

How to play with it

Hotkeys

1. Press SPACE to start or pause
2. “wasd” for moving the camera
3. When paused, press “p” to procceed only one step(useful for debugging)

Try different parameters

1. Try to change the parameters in “init_particles()” (like init_pos, cube size etc.) or the num_particles
2. Try to give another angle to the “rotation()”
3. Try to change the stiffness of penalty force in “collision_response()”

Switch the branches

The different branches record the progress of coding.
Start with the simplest case, and gradually add features

1. minimal_ggui: show two static particles in the screen
2. init_particles: neatly stacked particles to form a box
3. translation: move the box
4. rotation: rotate the box
5. collision_one_particle: colliding with a ramp, the collided particles are dyed to red
6. bounce: the box can drop and rebounce, but no rotation
7. master: the main branch

Add more features

1. Add more boundary walls (now there is only a ground)
2. Try to give better collision_reponse (more complex penalty forces or better)
3. Try to add collision detection for collision with complex geometries
4. Try to add more rigid bodies, and make them collide with each other

Want to understand the theory?

see:

• Matthias Müller, Bruno Heidelberger, Matthias Teschner, and Markus Gross. 2005. Meshless deformations based on shape matching. ACM Trans. Graph. 24, 3 (July 2005), 471–478. https://doi.org/10.1145/1073204.1073216
• GAMES103 Course, by Huamin Wang. https://games-cn.org/games103/

近期我会录个视频到B站从头讲解原理和实现方法。图文版也会放到博客里。

Codes

import taichi as ti
import math

ti.init()

num_particles = 2000
dim=3
world_scale_factor = 1.0/100.0
dt = 1e-2
mass_inv = 1.0

positions = ti.Vector.field(dim, float, num_particles)
velocities = ti.Vector.field(dim, float, num_particles)
pos_draw = ti.Vector.field(dim, float, num_particles)
force = ti.Vector.field(dim, float, num_particles)
penalty_force = ti.Vector.field(dim, float, num_particles)
positions0 = ti.Vector.field(dim, float, num_particles)
radius_vector = ti.Vector.field(dim, float, num_particles)
paused = ti.field(ti.i32, shape=())
q_inv = ti.Matrix.field(n=3, m=3, dtype=float, shape=())

@ti.kernel
def init_particles():
    init_pos = ti.Vector([70.0, 50.0, 0.0])
    cube_size = 20 
    spacing = 2 
    num_per_row = (int) (cube_size // spacing) + 1
    num_per_floor = num_per_row * num_per_row
    for i in range(num_particles):
        floor = i // (num_per_floor) 
        row = (i % num_per_floor) // num_per_row
        col = (i % num_per_floor) % num_per_row
        positions[i] = ti.Vector([col*spacing, floor*spacing, row*spacing]) + init_pos


@ti.kernel
def shape_matching():
    # update vel and pos firtly
    gravity = ti.Vector([0.0, -9.8, 0.0])
    for i in range(num_particles):
        positions0[i] = positions[i]
        f = gravity + penalty_force[i]
        velocities[i] += mass_inv * f * dt 
        positions[i] += velocities[i] * dt

    #compute the new(matched shape) mass center
    c = ti.Vector([0.0, 0.0, 0.0])
    for i in range(num_particles):
        c += positions[i]
    c /= num_particles

    #compute transformation matrix and extract rotation
    A = sum1 = ti.Matrix([[0.0] * 3 for _ in range(3)], ti.f32)
    for i in range(num_particles):
        sum1 += (positions[i] - c).outer_product(radius_vector[i])
    A = sum1 @ q_inv[None]

    R, _ = ti.polar_decompose(A)

    #update velocities and positions
    for i in range(num_particles):
        positions[i] = c + R @ radius_vector[i]
        velocities[i] = (positions[i] - positions0[i]) / dt


@ti.kernel
def compute_radius_vector():
    #compute the mass center and radius vector
    center_mass = ti.Vector([0.0, 0.0, 0.0])
    for i in range(num_particles):
        center_mass += positions[i]
    center_mass /= num_particles
    for i in range(num_particles):
        radius_vector[i] = positions[i] - center_mass


@ti.kernel
def precompute_q_inv():
    res = ti.Matrix([[0.0] * 3 for _ in range(3)], ti.f64)
    for i in range(num_particles):
        res += radius_vector[i].outer_product(radius_vector[i])
    q_inv[None] = res.inverse()


@ti.kernel
def collision_response():
    eps = 2.0 # the padding to prevent penatrating
    k = 20.0 # stiffness of the penalty force
    #boundary for skybox (xmin, ymin, zmin, xmax, ymax, zmax)
    boundary = ti.Matrix([[0.0, 0.0, 0.0], [100.0, 100.0, 100.0]], ti.f32)
    boundary[0,:] = boundary[0,:] + eps
    boundary[1,:] = boundary[1,:] - eps
    for i in range(num_particles):
        if positions[i].y < boundary[0,1]:
            n_dir = ti.Vector([0.0, 1.0, 0.0])
            phi = positions[i].y - boundary[0,1]
            penalty_force[i] = k * ti.abs(phi)  * n_dir
    
        #TODO: more boundaries...


@ti.kernel
def rotation(angle:ti.f32):
    theta = angle / 180.0 * math.pi
    R = ti.Matrix([
    [ti.cos(theta), -ti.sin(theta), 0.0], 
    [ti.sin(theta), ti.cos(theta), 0.0], 
    [0.0, 0.0, 1.0]
    ])
    for i in range(num_particles):
        positions[i] = [email protected][i]

# ---------------------------------------------------------------------------- #
# substep #
# ---------------------------------------------------------------------------- #
def substep():
    penalty_force.fill(0.0)
    collision_response()
    shape_matching()
# ---------------------------------------------------------------------------- #
# end substep #
# ---------------------------------------------------------------------------- #

@ti.kernel
def world_scale():
    for i in range(num_particles):
        pos_draw[i] = positions[i] * world_scale_factor

#init the window, canvas, scene and camerea
window = ti.ui.Window("rigidbody", (1024, 1024),vsync=True)
canvas = window.get_canvas()
scene = ti.ui.Scene()
camera = ti.ui.make_camera()

#initial camera position
camera.position(0.5, 1.0, 1.95)
camera.lookat(0.5, 0.3, 0.5)
camera.fov(55)

def main():
    init_particles()
    rotation(30)
    compute_radius_vector() #store the shape of rigid body
    precompute_q_inv()

    paused[None] = True
    while window.running:
        if window.get_event(ti.ui.PRESS):
            #press space to pause the simulation
            if window.event.key == ti.ui.SPACE:
                paused[None] = not paused[None]
            
            #proceed once for debug
            if window.event.key == 'p':
                substep()

        #do the simulation in each step
        if (paused[None] == False) :
            for i in range(int(0.05/dt)):
                substep()

        #set the camera, you can move around by pressing 'wasdeq'
        camera.track_user_inputs(window, movement_speed=0.03, hold_key=ti.ui.RMB)
        scene.set_camera(camera)

        #set the light
        scene.point_light(pos=(0, 1, 2), color=(1, 1, 1))
        scene.point_light(pos=(0.5, 1.5, 0.5), color=(0.5, 0.5, 0.5))
        scene.ambient_light((0.5, 0.5, 0.5))
        
        #draw particles
        world_scale()
        scene.particles(pos_draw, radius=0.01, color=(0, 1, 1))

        #show the frame
        canvas.scene(scene)
        window.show()

if __name__ == '__main__':
    main()