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

Github repo

https://github.com/chunleili/tiRigidBody

Demo

在这里插入图片描述

原理与代码讲解

视频讲解

https://www.bilibili.com/video/BV1tr4y1L7qi/)

图文讲解

(93条消息) 【物理模拟】最简单的shape matching的原理与实践_beidou111的博客-CSDN博客

How to run

First, intall latest taichi by

pip install taichi

Then, run with

python rigidbody.py

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:

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] = R@positions[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()

7 个赞

更新一下。之前地板碰撞用的罚函数,貌似有点小BUG。下面给改正了。改动部分只有几行。

这是改动后的效果

demo

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)
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
        velocities[i] += mass_inv * f * dt 
        positions[i] += velocities[i] * dt
        if positions[i].y < 0.0:
            positions[i] = positions0[i]
            positions[i].y = 0.0
 
    #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 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] = R@positions[i]

# ---------------------------------------------------------------------------- #
#                                    substep                                   #
# ---------------------------------------------------------------------------- #
def substep():
    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", (1280, 720),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(60) #initially rotate the cube
    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()
3 个赞

请问一下 你这个旋转轴是怎么选取呢?

老师您好:
想请问一下,如果我想模拟任意形状的刚体运动,我该如何生成规则网格数据?我目前没有找到有关于二维形状规则网格的划分的工具,请问有什么好的办法吗?非常感谢!