课件与录像汇总

最新网页请点击这里 (2020年11月12日,不断更新)

Lec 1 导论

Lec 2 拉格朗日视角(1)

Lec 3 拉格朗日视角(2)

Lec 4 欧拉视角

Lec 5 多体问题与涡方法

Lec 6 有限元与拓扑优化

Lec 7 混合欧拉-拉格朗日视角 (1)

Lec 8 混合欧拉-拉格朗日视角 (2)

Lec 9 高性能物理模拟

Lec 10 大结局

18 个赞

原谅我太爱fractal了,看到就忍不住了。 :smiling_face_with_three_hearts:谢谢胡老师送的六一礼物

import taichi as ti
ti.init(arch=ti.gpu)
n = 320
pixels = ti.var(dt=ti.f32, shape=(n * 2, n))

@ti.func
def complex_power(z, power: ti.i32):
    r = ti.sqrt(z[0]**2 + z[1]**2)
    theta = ti.atan2(z[1], z[0])
    return ti.Vector([r**power * ti.cos(power*theta), r**power * ti.sin(power*theta)])

@ti.kernel
def paint(t: ti.f32, power: ti.i32):
    for i, j in pixels:  # Parallized over all pixels
        # Julia Set
        freq = 1.0 / power
        c = ti.Vector([0.7885 * ti.cos(freq*t), 0.7885 * ti.sin(freq*t)])
        z = ti.Vector([i / n - 1, j / n - 0.5]) * 2

        iterations = 0
        while z.norm() < 20 and iterations < 50:
            z = complex_power(z, power) + c
            iterations += 1
        pixels[i, j] = 1 - iterations * 0.02

power = eval(input("Power of z -> "))
gui = ti.GUI("Julia Set", res=(n * 2, n))

for i in range(1000000):
    paint(i * 0.03, power)
    gui.set_image(pixels)
    gui.show()
8 个赞

有意思呀!我跑了你的程序z=1,2,3,4,5,6,7,...,原来每个数看起来结果都不太一样。可否单独开一个帖子,顺便贴几个图?

1 个赞

是的,fractal根据z的power不同出现在实数部分的样子不同。 明天看看有空的话写个帖子,今晚先睡了 :partying_face:

3 个赞

-1 -2 -3 更有意思!完全不是传统的分形图案了,反而是一堆小球在运动

3 个赞

太神奇了

请问胡老师会讲解 Giga-Voxel Narrowband TopOpt 代码吗,谢谢!

lec3 的视频呢哭唧唧

后面应该会讲一讲简单的linear elasticity FEM (quad mesh),也会提到multigrid和topology optimization,但是那篇文章里面的算法有点硬核,所以不会过于深入。

Lec3的课件lec3-elasticity中Neo-Hookean的P(F)错了,应该是 F - F-T
附一下基于demo FEM99的手动求导版本 ,感谢胡老师超棒的课程和超棒的taichi : )

import taichi as ti
ti.init(arch=ti.gpu)

N = 16
dt = 1e-4
dx = 1 / N
rho = 4e1
NF = 2 * N ** 2   # number of faces
NV = (N + 1) ** 2 # number of vertices
E, nu = 4e4, 0.2  # Young's modulus and Poisson's ratio
mu, lam = E / 2 / (1 + nu), E * nu / (1 + nu) / (1 - 2 * nu) # Lame parameters
ball_pos, ball_radius = ti.Vector([0.5, 0.0]), 0.32
gravity = ti.Vector([0, -40])
damping = 12.5

pos = ti.Vector.field(2, float, NV)
vel = ti.Vector.field(2, float, NV)
f2v = ti.Vector.field(3, int, NF)  # ids of three vertices of each face
B = ti.Matrix.field(2, 2, float, NF)
W = ti.field(float, NF)
phi = ti.field(float, NF)  # potential energy of each face (Neo-Hookean)
U = ti.field(float, (), needs_grad=True)  # total potential energy
f = ti.Vector.field(2, float, NV)

@ti.kernel
def update_force():
    for i in range(NF):
        ia, ib, ic = f2v[i]
        a, b, c = pos[ia], pos[ib], pos[ic]

        D_i = ti.Matrix.cols([a - c, b - c])
        F = D_i @ B[i]
        F_it =  F.inverse().transpose()

        PF = mu*(F - F_it) + lam* ti.log(F.determinant())*F_it
        H = -W[i]* PF@B[i].transpose()

        fa = ti.Vector([H[0, 0], H[1, 0]])
        fb = ti.Vector([H[0, 1], H[1, 1]])
        fc = -fa-fb
        f[ia] += fa
        f[ib] += fb
        f[ic] += fc

@ti.kernel
def advance():
    for i in range(NV):
        acc = f[i] / (rho * dx ** 2)
        vel[i] += dt * (acc + gravity)
        vel[i] *= ti.exp(-dt * damping)
    for i in range(NV):
        # ball boundary condition:
        disp = pos[i] - ball_pos
        disp2 = disp.norm_sqr()
        if disp2 <= ball_radius ** 2:
            NoV = vel[i].dot(disp)
            if NoV < 0: vel[i] -= NoV * disp / disp2
        # rect boundary condition:
        cond = pos[i] < 0 and vel[i] < 0 or pos[i] > 1 and vel[i] > 0
        for j in ti.static(range(pos.n)):
           if cond[j]: vel[i][j] = 0
        pos[i] += dt * vel[i]

@ti.kernel
def init_pos():
    for i, j in ti.ndrange(N + 1, N + 1):
        k = i * (N + 1) + j
        pos[k] = ti.Vector([i, j]) / N * 0.25 + ti.Vector([0.45, 0.45])
        vel[k] = ti.Vector([0, 0])
    for i in range(NF):
        ia, ib, ic = f2v[i]
        a, b, c = pos[ia], pos[ib], pos[ic]
        B_i = ti.Matrix.cols([a - c, b - c])
        B[i] = B_i.inverse()
        W[i] = abs( B_i.determinant())

@ti.kernel
def init_mesh():
    for i, j in ti.ndrange(N, N):
        k = (i * N + j) * 2
        a = i * (N + 1) + j
        b = a + 1
        c = a + N + 2
        d = a + N + 1
        f2v[k + 0] = [a, b, c]
        f2v[k + 1] = [c, d, a]

@ti.kernel
def init_parameter():
    for i in f:
        f[i] = ti.Vector([0.0, 0.0])

init_mesh()
init_pos()
gui = ti.GUI('FEM99')
while gui.running:
    for e in gui.get_events():
        if e.key == gui.ESCAPE:
            gui.running = False
        elif e.key == 'r':
            init_pos()
    for i in range(30):
        init_parameter()
        update_force()
        advance()
    gui.circles(pos.to_numpy(), radius=2, color=0xffaa33)
    gui.circle(ball_pos, radius=ball_radius * 512, color=0x666666)
    gui.show()
4 个赞