How to use the taichi.Matrix function in taichi
To help you get started, we’ve selected a few taichi examples, based on popular ways it is used in public projects.
Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
yuanming-hu / taichi / tests / python / test_dynamic.py View on Github 
def test_dynamic_matrix():
return
if ti.cfg.arch == ti.cuda:
return
x = ti.Matrix(3, 2, dt=ti.f32)
n = 8192
@ti.layout
def place():
ti.root.dynamic(ti.i, n, chunk_size=128).place(x)
@ti.kernel
def func():
ti.serialize()
for i in range(n // 4):
x[i * 4][1, 0] = i
func()
for i in range(n // 4):
assert x[i * 4][1, 0] == iself.inv_dx = 1.0 / self.dx
self.default_dt = 2e-2 * self.dx / size * dt_scale
self.p_vol = self.dx**self.dim
self.p_rho = 1000
self.p_mass = self.p_vol * self.p_rho
self.max_num_particles = max_num_particles
self.gravity = ti.Vector(self.dim, dt=ti.f32, shape=())
self.source_bound = ti.Vector(self.dim, dt=ti.f32, shape=2)
self.source_velocity = ti.Vector(self.dim, dt=ti.f32, shape=())
self.pid = ti.var(ti.i32)
# position
self.x = ti.Vector(self.dim, dt=ti.f32)
# velocity
self.v = ti.Vector(self.dim, dt=ti.f32)
# affine velocity field
self.C = ti.Matrix(self.dim, self.dim, dt=ti.f32)
# deformation gradient
self.F = ti.Matrix(self.dim, self.dim, dt=ti.f32)
# material id
self.material = ti.var(dt=ti.i32)
self.color = ti.var(dt=ti.i32)
# plastic deformation volume ratio
self.Jp = ti.var(dt=ti.f32)
if self.dim == 2:
indices = ti.ij
else:
indices = ti.ijk
offset = tuple(-self.grid_size // 2 for _ in range(self.dim))
self.offset = offsetdef _test_polar_decomp(dim, dt):
m = ti.Matrix(dim, dim, dt)
r = ti.Matrix(dim, dim, dt)
s = ti.Matrix(dim, dim, dt)
I = ti.Matrix(dim, dim, dt)
D = ti.Matrix(dim, dim, dt)
@ti.layout
def place():
ti.root.place(m, r, s, I, D)
@ti.kernel
def polar():
R, S = ti.polar_decompose(m[None], dt)
r[None] = R
s[None] = S
m[None] = R @ S
I[None] = R @ ti.transposed(R)
D[None] = S - ti.transposed(S)
def V(i, j):
return i * 2 + j * 7 + int(i == j) * 3def 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=dim, val=1) + dt * C[f, p]) @ F[f, p]
J = ti.determinant(new_F)
if particle_type[p] == 0: # fluid
sqrtJ = ti.sqrt(J)
# TODO: need pow(x, 1/3)
new_F = ti.Matrix([[sqrtJ, 0, 0], [0, sqrtJ, 0], [0, 0, 1]])
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
# ti.print(act)
A = ti.Matrix([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) * act
cauchy = ti.Matrix(zero_matrix())
mass = 0.0
ident = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
if particle_type[p] == 0: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
# ti.print(act)
A = ti.Matrix([[0.0, 0.0], [0.0, 1.0]]) * act
cauchy = ti.Matrix([[0.0, 0.0], [0.0, 0.0]])
mass = 0.0
if particle_type[p] == 0:
mass = 4
cauchy = ti.Matrix([[1.0, 0.0], [0.0, 0.1]]) * (J - 1) * E
else:
mass = 1
cauchy = 2 * mu * (new_F - r) @ ti.transposed(new_F) + \
ti.Matrix.diag(2, la * (J - 1) * J)
cauchy += new_F @ A @ ti.transposed(new_F)
stress = -(dt * p_vol * 4 * inv_dx * inv_dx) * cauchy
affine = stress + 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].atomic_add(
weight * (mass * v[f, p] + affine @ dpos))
grid_m_in[base + offset].atomic_add(weight * mass)act = actuation[f, ti.max(0, act_id)] * act_strength
if act_id == -1:
act = 0.0
# ti.print(act)
A = ti.Matrix([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]]) * act
cauchy = ti.Matrix(zero_matrix())
mass = 0.0
ident = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
if particle_type[p] == 0:
mass = 4
cauchy = ti.Matrix(ident) * (J - 1) * E
else:
mass = 1
cauchy = mu * (new_F @ ti.transposed(new_F)) + ti.Matrix(ident) * (
la * ti.log(J) - mu)
cauchy += new_F @ A @ ti.transposed(new_F)
stress = -(dt * p_vol * 4 * inv_dx * inv_dx) * cauchy
affine = stress + mass * C[f, p]
for i in ti.static(range(3)):
for j in ti.static(range(3)):
for k in ti.static(range(3)):
offset = ti.Vector([i, j, k])
dpos = (ti.cast(ti.Vector([i, j, k]), real) - fx) * dx
weight = w[i](0) * w[j](1) * w[k](2)
grid_v_in[base + offset].atomic_add(
weight * (mass * v[f, p] + affine @ dpos))
grid_m_in[base + offset].atomic_add(weight * mass)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
# ti.print(act)
A = ti.Matrix([[0.0, 0.0], [0.0, 1.0]]) * act
cauchy = ti.Matrix([[0.0, 0.0], [0.0, 0.0]])
mass = 0.0
if particle_type[p] == 0:
mass = 4
cauchy = ti.Matrix([[1.0, 0.0], [0.0, 0.1]]) * (J - 1) * E
else:
mass = 1
cauchy = 2 * mu * (new_F - r) @ ti.transposed(new_F) + \
ti.Matrix.diag(2, la * (J - 1) * J)
cauchy += new_F @ A @ ti.transposed(new_F)
stress = -(dt * p_vol * 4 * inv_dx * inv_dx) * cauchy
affine = stress + 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].atomic_add(
weight * (mass * v[f, p] + affine @ dpos))
grid_m_in[base + offset].atomic_add(weight * mass)def rotation_matrix(r):
return ti.Matrix([[ti.cos(r), -ti.sin(r)], [ti.sin(r), ti.cos(r)]])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
