Homework2: 整合mpm_lagrangian_forces和mpm_99

其实就是把 mpm_99 复制粘贴到了 mpm_lagrangian_forces 里面,抄了一些导数的式子,免去了自动求导。
去掉自动求导之后似乎快了不少。
实验结果见 https://www.bilibili.com/video/BV1SV411z76R/

实现了三种边界情况,可以看到 sticky 的边界条件下水没法沿着边界向上滑动,靠近边界的物质看起来也不怎么容易移动。
separate 与 slip 边界条件很相像,但是因为 slip 边界条件不允许在边界上出现法向量方向的速度,水珠在 y 的上边界(或者说天花板)不会怎么受到重力作用,下落速度大大减慢,表现出在天花板上滑动。

References:
GAMES201
mpm_99, mpm_lagrangian_forces (from ti example)
Barbic and Sifakis’ SIGGRAPH Course notes on FEM / Model Reduction (2012)
Jiang et al. SIGGRAPH Course: The Material Point Method for Simulating Continuum Materials (2016)

import taichi as ti, sys
ti.init(debug = False, arch = ti.gpu)

bc_mode = int(sys.argv[1]) if len(sys.argv) > 1 else 2 # 0 for sticky; 1 for slip; 2 for separate

real = ti.f32
dim = 2
n_particle_x = 60
n_particle_y = 6
n_each = 2048
n_material = 3
n_mesh = n_particle_x * n_particle_y
n_particles = n_mesh + n_each * n_material
n_elements = (n_particle_x - 1) * (n_particle_y - 1) * 2
n_grid = 128
dx = 1 / n_grid
dx_inv = 1 / dx
dt = 2e-5
n_subframes = int(2e-3 // dt)
rho_p = 1
vol_p = (dx * 0.5) ** 2
m_p = vol_p * rho_p
E, nu = 5e3, 0.2
mu_0, lambda_0 = E / (2 * (1 + nu)), E * nu / ((1 + nu) * (1 - 2 * nu))  # Lame parameters
x_mesh_0, y_mesh_0 = 0.3, 0.9
dx_mesh = dx
bound = 3

x = ti.Vector.field(2, dtype = ti.f32, shape = n_particles)
v = ti.Vector.field(2, dtype = ti.f32, shape = n_particles)
f = ti.Vector.field(2, dtype = ti.f32, shape = n_particles)
C = ti.Matrix.field(2, 2, dtype = ti.f32, shape = n_particles)
Jp = ti.field(dtype = ti.f32, shape = n_particles)
material = ti.field(dtype = ti.i32, shape = n_particles)
vertices = ti.Vector.field(3, dtype = ti.i32, shape = n_elements)
B = ti.Matrix.field(2, 2, dtype = ti.f32, shape = n_elements)
WD = ti.Matrix.field(2, 2, dtype = ti.f32, shape = n_elements)
F = ti.Matrix.field(2, 2, dtype = ti.f32, shape = n_particles)
v_grid = ti.Vector.field(2, dtype = ti.f32, shape = (n_grid, n_grid))
m_grid = ti.field(dtype = ti.f32, shape = (n_grid, n_grid))
gravity = ti.Vector.field(2, dtype = ti.f32, shape = ())

@ti.func
def ij2p(i, j):
    return i * ti.static(n_particle_y) + j

@ti.func
def get_T(vertices):
    a = x[vertices[0]] - x[vertices[2]]
    b = x[vertices[1]] - x[vertices[2]]
    return ti.Matrix([[a.x, b.x], [a.y, b.y]])

@ti.kernel
def init():
    gravity[None] = ti.Vector([0, -9.8 * 10])
    for i, j in ti.ndrange(n_particle_x, n_particle_y):
        p = ij2p(i, j)
        x[p] = [x_mesh_0 + i * dx_mesh, y_mesh_0 + j * dx_mesh]
        v[p] = ti.Vector.zero(ti.f32, 2)
        C[p] = ti.Matrix.identity(ti.f32, 2)
        material[p] = 3
    for i, j in ti.ndrange(n_particle_x - 1, n_particle_y - 1):
        e = i * (n_particle_y - 1) + j
        vertices[e * 2] = ti.Vector([ij2p(i, j), ij2p(i + 1, j), ij2p(i, j + 1)])
        tmp = get_T(vertices[e * 2])
        B[e * 2] = tmp.inverse()
        WD[e * 2] = -ti.abs(tmp.determinant()) / 2 * tmp.inverse().transpose()
        vertices[e * 2 + 1] = ti.Vector([ij2p(i + 1, j + 1), ij2p(i, j + 1), ij2p(i + 1, j)])
        tmp = get_T(vertices[e * 2 + 1])
        B[e * 2 + 1] = tmp.inverse()
        WD[e * 2 + 1] = -ti.abs(tmp.determinant()) / 2 * tmp.inverse().transpose()
    for p in range(n_each * n_material):
        x[p + n_mesh] = [ti.random() * 0.2 + 0.1 + 0.1 * (p // n_each), ti.random() * 0.2 + 0.1 + 0.21 * (p // n_each)]
        v[p + n_mesh] = [0, 0]
        F[p + n_mesh] = ti.Matrix.identity(ti.f32, 2)
        Jp[p + n_mesh] = 1
        C[p + n_mesh] = ti.Matrix.zero(ti.f32, 2, 2)
        material[p + n_mesh] = p // n_each

@ti.kernel
def get_force():
    for e in vertices:
        F_tmp = get_T(vertices[e]) @ B[e]
        J = F_tmp.determinant()
        mu, la = mu_0 * 1e-3, lambda_0 * 1e-3
        P = mu * (F_tmp - F_tmp.transpose().inverse()) + la * ti.log(J) * F_tmp.transpose().inverse()
        H = P @ WD[e] * 1e-3 * (8 * n_grid ** 2) # TODO: magic number
        tmp_x = ti.Vector([H[0, 0], H[1, 0]])
        tmp_y = ti.Vector([H[0, 1], H[1, 1]])
        f[vertices[e].x] += tmp_x
        f[vertices[e].y] += tmp_y
        f[vertices[e].z] += -tmp_x - tmp_y

@ti.kernel
def p2g():
    for p in x:
        base = (x[p] * dx_inv - 0.5).cast(int)
        local = x[p] * dx_inv - base
        w = [0.5 * (1.5 - local) ** 2, 0.75 - (local - 1) ** 2, 0.5 * (local - 0.5) ** 2]
        affine = m_p * C[p]
        if material[p] != 3:
            F[p] = (ti.Matrix.identity(ti.f32, 2) + dt * C[p]) @ F[p]
            U, sig, V = ti.svd(F[p])
            if material[p] == 0:
                J = sig[0, 0] * sig[1, 1]
                F[p] = ti.Matrix.identity(ti.f32, 2) * ti.sqrt(J)
                affine += -4 * dt * (dx_inv ** 2) * vol_p * E * (J - 1) * ti.Matrix.identity(ti.f32, 2)
            elif material[p] == 1:
                J = sig[0, 0] * sig[1, 1]
                h = 0.3
                mu, la = mu_0 * h, lambda_0 * h
                stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose() + la * (J - 1) * J * ti.Matrix.identity(ti.f32, 2)
                affine += -4 * dt * (dx_inv ** 2) * vol_p * stress
            elif material[p] == 2:
                h = max(0.1, min(5, ti.exp(10 * (1.0 - Jp[p]))))  # Hardening coefficient: snow gets harder when compressed
                mu, la = mu_0 * h, lambda_0 * h
                J = 1.0
                for d in ti.static(range(2)):
                    new_sig = min(max(sig[d, d], 1 - 2.5e-2), 1 + 4.5e-3)  # Plasticity
                    Jp[p] *= sig[d, d] / new_sig
                    sig[d, d] = new_sig
                    J *= new_sig
                F[p] = U @ sig @ V.transpose()
                stress = 2 * mu * (F[p] - U @ V.transpose()) @ F[p].transpose() + la * (J - 1) * J * ti.Matrix.identity(ti.f32, 2)
                affine += -4 * dt * (dx_inv ** 2) * vol_p * stress
        for i, j in ti.static(ti.ndrange(3, 3)):
            offset = ti.Vector([i, j])
            k = w[i].x * w[j].y
            dpos = (offset.cast(float) - local) * dx
            v_grid[base + offset] += k * (m_p * v[p] + affine @ dpos + m_p * f[p])
            m_grid[base + offset] += k * m_p

@ti.kernel
def grid_op():
    for i, j in v_grid:
        I = ti.Vector([i, j])
        if m_grid[I] <= 0:
            continue
        v_grid[I] = v_grid[I] / m_grid[I] + gravity * dt
        if ti.static(bc_mode) == 0: # sticky
            if i < bound or i >= n_grid - bound or j < bound or j >= n_grid - bound: v_grid[I] = [0, 0]
        elif ti.static(bc_mode) == 1: # slip
            if i < bound or i >= n_grid - bound: v_grid[I].x = 0
            if j < bound or j >= n_grid - bound: v_grid[I].y = 0
        elif ti.static(bc_mode) == 2: # separate
            if i < bound and v_grid[I].x < 0: v_grid[I].x = 0
            if i >= n_grid - bound and v_grid[I].x > 0: v_grid[I].x = 0
            if j < bound and v_grid[I].y < 0: v_grid[I].y = 0
            if j >= n_grid - bound and v_grid[I].y > 0: v_grid[I].y = 0

@ti.kernel
def g2p():
    for p in x:
        base = (x[p] * dx_inv - 0.5).cast(int)
        local = x[p] * dx_inv - base
        w = [0.5 * (1.5 - local) ** 2, 0.75 - (local - 1) ** 2, 0.5 * (local - 0.5) ** 2]
        new_v = ti.Vector.zero(ti.f32, 2)
        new_C = ti.Matrix.zero(ti.f32, 2, 2)
        for i, j in ti.static(ti.ndrange(3, 3)):
            k = w[i].x * w[j].y
            offset = ti.Vector([i, j])
            dpos = offset.cast(float) - local # actually dpos / dx
            new_v += k * v_grid[base + offset]
            new_C += k * 4 * dx_inv * v_grid[base + offset].outer_product(dpos)
        x[p] += dt * new_v
        v[p] = new_v
        C[p] = new_C

def substep():
    v_grid.fill(0)
    m_grid.fill(0)
    f.fill(0)
    get_force()
    p2g()
    grid_op()
    g2p()

def main():
    gui = ti.GUI('trivial mpm with lagrangian force', res = (512, 512), background_color = 0x112F41)
    result_dir = "./results"
    #video_manager = ti.VideoManager(output_dir=result_dir, framerate=24, automatic_build=False)
    init()
    breakflag = False
    while gui.running:
        for e in gui.get_events(ti.GUI.PRESS):
            if e.key in [ti.GUI.ESCAPE, ti.GUI.EXIT]: breakflag = True
            elif e.key == 'r': init()
            elif e.key == 'g':
                for subframe in range(n_subframes):
                    substep()
        if breakflag: break
        for subframe in range(n_subframes):
            substep()
        gui.circles(x.to_numpy()[0: n_mesh], radius = 1, color = 0xF2B134)
        gui.circles(x.to_numpy()[n_mesh: n_mesh + n_each], color = 0x068587)
        gui.circles(x.to_numpy()[n_mesh + n_each: n_mesh + n_each * 2], color = 0xED553B)
        gui.circles(x.to_numpy()[n_mesh + n_each * 2: n_particles], color = 0xEEEEF0)
        #video_manager.write_frame(gui.get_image())
        gui.show()
    #video_manager.make_video(gif = True, mp4 = True)

if __name__ == '__main__': main()

3 Likes

大赞!!! :100: 做了个GIF方便大家分享:

video

1 Like