@ti.kernel
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=dim, val=1) + (dt * C[(f, p)])) @ F[(f, p)])
        J = ti.determinant(new_F)
        if (particle_type[p] == 0):
            sqrtJ = ti.sqrt(J)
            new_F = ti.Matrix([[sqrtJ, 0, 0], [0, sqrtJ, 0], [0, 0, 1]])
        F[((f + 1), p)] = new_F
        act_id = actuator_id[p]
        act = (actuation[(f, ti.max(0, act_id))] * act_strength)
        if (act_id == (- 1)):
            act = 0.0
        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))
@ti.kernel
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=dim, val=1) + (dt * C[(f, p)])) @ F[(f, p)])
        J = ti.determinant(new_F)
        if (particle_type[p] == 0):
            sqrtJ = ti.sqrt(J)
            new_F = ti.Matrix([[sqrtJ, 0, 0], [0, sqrtJ, 0], [0, 0, 1]])
        F[((f + 1), p)] = new_F
        act_id = actuator_id[p]
        act = (actuation[(f, ti.max(0, act_id))] * act_strength)
        if (act_id == (- 1)):
            act = 0.0
        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))