Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def g2p(f: ti.i32):
for p in range(0, n_particles):
base = ti.cast(x[f, p] * inv_dx - 0.5, ti.i32)
fx = x[f, p] * inv_dx - ti.cast(base, real)
w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1.0),
0.5 * ti.sqr(fx - 0.5)]
new_v = ti.Vector([0.0, 0.0])
new_C = ti.Matrix([[0.0, 0.0], [0.0, 0.0]])
for i in ti.static(range(3)):
for j in ti.static(range(3)):
dpos = ti.cast(ti.Vector([i, j]), real) - fx
g_v = grid_v_out[base(0) + i, base(1) + j]
weight = w[i](0) * w[j](1)
new_v += weight * g_v
new_C += 4 * weight * ti.outer_product(g_v, dpos) * inv_dx
v[f + 1, p] = new_v
x[f + 1, p] = x[f, p] + dt * v[f + 1, p]
C[f + 1, p] = new_C
def render_photon_map(t: ti.i32, offset_x: ti.f32, offset_y: ti.f32):
for i in range(n_grid): # Parallelized over GPU threads
for j in range(n_grid):
grad = height_gradient[i, j] * (1 - offset_x) * (1 - offset_y) + \
height_gradient[i + 1, j] * offset_x * (1 - offset_y) + \
height_gradient[i, j + 1] * (1 - offset_x) * offset_y + \
height_gradient[i + 1, j + 1] * offset_x * offset_y
scale = 5.0
sample_x = i - grad[0] * scale + offset_x
sample_y = j - grad[1] * scale + offset_y
sample_x = ti.min(n_grid - 1, ti.max(0, sample_x))
sample_y = ti.min(n_grid - 1, ti.max(0, sample_y))
sample_xi = ti.cast(ti.floor(sample_x), ti.i32)
sample_yi = ti.cast(ti.floor(sample_y), ti.i32)
frac_x = sample_x - sample_xi
frac_y = sample_y - sample_yi
x = sample_xi
y = sample_yi
ti.atomic_add(rendered[x, y], (1 - frac_x) * (1 - frac_y))
ti.atomic_add(rendered[x, y + 1], (1 - frac_x) * frac_y)
ti.atomic_add(rendered[x + 1, y], frac_x * (1 - frac_y))
ti.atomic_add(rendered[x + 1, y + 1], frac_x * frac_y)
rotated_x = dir[0] * ti.cos(angle) + dir[2] * ti.sin(angle)
rotated_z = -dir[0] * ti.sin(angle) + dir[2] * ti.cos(angle)
dir[0] = rotated_x
dir[2] = rotated_z
point = camera_origin + (k + 1) * dx * dir
# Convert to coordinates of the density grid box
box_x = point[0] + 0.5
box_y = point[1] + 0.5
box_z = point[2] + 0.5
# Density grid location
index_x = ti.cast(ti.floor(box_x * density_res), ti.i32)
index_y = ti.cast(ti.floor(box_y * density_res), ti.i32)
index_z = ti.cast(ti.floor(box_z * density_res), ti.i32)
index_x = ti.max(0, ti.min(index_x, density_res - 1))
index_y = ti.max(0, ti.min(index_y, density_res - 1))
index_z = ti.max(0, ti.min(index_z, density_res - 1))
flag = 0
if in_box(point[0], point[1], point[2]):
flag = 1
contribution = density[index_z, index_y, index_x] * flag
ti.atomic_add(field[view_id, y, x], contribution)
def g2p(f: ti.i32):
for p in range(n_particles):
base = ti.cast(x[f, p] * inv_dx - 0.5, ti.i32)
fx = x[f, p] * inv_dx - ti.cast(base, real)
w = [
0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1.0), 0.5 * ti.sqr(fx - 0.5)
]
new_v = ti.Vector([0.0, 0.0])
new_C = ti.Matrix([[0.0, 0.0], [0.0, 0.0]])
for i in ti.static(range(3)):
for j in ti.static(range(3)):
dpos = ti.cast(ti.Vector([i, j]), real) - fx
g_v = grid_v_out[f, base(0) + i, base(1) + j]
weight = w[i](0) * w[j](1)
new_v += weight * g_v
new_C += 4 * weight * ti.outer_product(g_v, dpos) * inv_dx
v[f + 1, p] = new_v
x[f + 1, p] = x[f, p] + dt * v[f + 1, p]
def p2g():
for p in x:
base = ti.cast(x[p] * inv_dx - 0.5, ti.i32)
fx = x[p] * inv_dx - ti.cast(base, ti.f32)
w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1), 0.5 * ti.sqr(fx - 0.5)]
affine = p_mass * C[p]
for i in ti.static(range(3)):
for j in ti.static(range(3)):
offset = ti.Vector([i, j])
dpos = (ti.cast(ti.Vector([i, j]), ti.f32) - fx) * dx
weight = w[i](0) * w[j](1)
grid_v[base + offset].atomic_add(
weight * (p_mass * v[p] - x.grad[p] + affine @ dpos))
grid_m[base + offset].atomic_add(weight * p_mass)
def p2g(f: ti.i32):
for p in range(0, n_particles):
base = ti.cast(x[f, p] * inv_dx - 0.5, ti.i32)
fx = x[f, p] * inv_dx - ti.cast(base, ti.i32)
w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1), 0.5 * ti.sqr(fx - 0.5)]
new_F = (ti.Matrix.diag(dim=2, val=1) + dt * C[f, p]) @ F[f, p]
F[f + 1, p] = new_F
J = ti.determinant(new_F)
r, s = ti.polar_decompose(new_F)
cauchy = 2 * mu * (new_F - r) @ ti.transposed(new_F) + \
ti.Matrix.diag(2, la * (J - 1) * J)
stress = -(dt * p_vol * 4 * inv_dx * inv_dx) * cauchy
affine = stress + p_mass * C[f, p]
for i in ti.static(range(3)):
for j in ti.static(range(3)):
offset = ti.Vector([i, j])
dpos = (ti.cast(ti.Vector([i, j]), real) - fx) * dx
weight = w[i](0) * w[j](1)
grid_v_in[f, base + offset] += weight * (
def p2g(f: ti.i32):
for p in range(0, n_particles):
base = ti.cast(x[f, p] * inv_dx - 0.5, ti.i32)
fx = x[f, p] * inv_dx - ti.cast(base, ti.i32)
w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1), 0.5 * ti.sqr(fx - 0.5)]
new_F = (ti.Matrix.diag(dim=2, val=1) + dt * C[f, p]) @ F[f, p]
F[f + 1, p] = new_F
J = ti.determinant(new_F)
r, s = ti.polar_decompose(new_F)
cauchy = 2 * mu * (new_F - r) @ ti.transposed(new_F) + \
ti.Matrix.diag(2, la * (J - 1) * J)
stress = -(dt * p_vol * 4 * inv_dx * inv_dx) * cauchy
affine = stress + p_mass * C[f, p]
for i in ti.static(range(3)):
for j in ti.static(range(3)):
offset = ti.Vector([i, j])
dpos = (ti.cast(ti.Vector([i, j]), real) - fx) * dx
weight = w[i](0) * w[j](1)
grid_v_in[base + offset] += weight * (p_mass * v[f, p] + affine @ dpos)
def p2g(f: ti.i32):
for p in range(0, n_particles):
base = ti.cast(x[f, p] * inv_dx - 0.5, ti.i32)
fx = x[f, p] * inv_dx - ti.cast(base, ti.i32)
w = [0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1), 0.5 * ti.sqr(fx - 0.5)]
new_F = (ti.Matrix.diag(dim=2, val=1) + dt * C[f, p]) @ F[f, p]
J = ti.determinant(new_F)
if particle_type[p] == 0: # fluid
sqrtJ = ti.sqrt(J)
new_F = ti.Matrix([[sqrtJ, 0], [0, sqrtJ]])
F[f + 1, p] = new_F
r, s = ti.polar_decompose(new_F)
act_id = actuator_id[p]
act = actuation[f, ti.max(0, act_id)] * act_strength
if act_id == -1:
act = 0.0
def g2p(f: ti.i32):
for p in range(0, n_particles):
base = ti.cast(x[f, p] * inv_dx - 0.5, ti.i32)
fx = x[f, p] * inv_dx - ti.cast(base, real)
w = [
0.5 * ti.sqr(1.5 - fx), 0.75 - ti.sqr(fx - 1.0), 0.5 * ti.sqr(fx - 0.5)
]
new_v = ti.Vector(zero_vec())
new_C = ti.Matrix(zero_matrix())
for i in ti.static(range(3)):
for j in ti.static(range(3)):
for k in ti.static(range(3)):
dpos = ti.cast(ti.Vector([i, j, k]), real) - fx
g_v = grid_v_out[base(0) + i, base(1) + j, base(2) + k]
weight = w[i](0) * w[j](1) * w[k](2)
new_v += weight * g_v
new_C += 4 * weight * ti.outer_product(g_v, dpos) * inv_dx
v[f + 1, p] = new_v
x[f + 1, p] = x[f, p] + dt * v[f + 1, p]
C[f + 1, p] = new_C
def advect(field: ti.template(), field_out: ti.template(),
t_offset: ti.template(), t: ti.i32):
"""Move field smoke according to x and y velocities (vx and vy)
using an implicit Euler integrator."""
for y in range(n_grid):
for x in range(n_grid):
center_x = y - v[t + t_offset, y, x][0]
center_y = x - v[t + t_offset, y, x][1]
# Compute indices of source cell
left_ix = ti.cast(ti.floor(center_x), ti.i32)
top_ix = ti.cast(ti.floor(center_y), ti.i32)
rw = center_x - left_ix # Relative weight of right-hand cell
bw = center_y - top_ix # Relative weight of bottom cell
# Wrap around edges
# TODO: implement mod (%) operator
left_ix = imod(left_ix, n_grid)
right_ix = inc_index(left_ix)
top_ix = imod(top_ix, n_grid)
bot_ix = inc_index(top_ix)
# Linearly-weighted sum of the 4 surrounding cells
field_out[t, y, x] = (1 - rw) * (
(1 - bw) * field[t - 1, left_ix, top_ix] +
bw * field[t - 1, left_ix, bot_ix]) + rw * (
(1 - bw) * field[t - 1, right_ix, top_ix] +