Hello!我写了一个简单的函数,输入四个Vector,然后输出一个Matrix。然而在给Matrix元素赋值的步骤报错如下。代码附在了后面。似乎是40行这一步,k[2 * i, 2 * j]
的index有一些问题。请问应该如何修改?谢谢!
另外,我这个函数k_ele()
目的是在kernel里结构for循环中建立一个8*8矩阵。请问有没有什么更好的方法实现?谢谢!
代码:
import taichi as ti
import numpy as np
ti.init(arch = ti.gpu)
# material parameters
nu = 0.3 #Poisson's ratio
l = 8.0 #length
# simulation parameters
n_no_x = 3 #no. of nodes in x direction
n_no_y = 2 #no. of nodes in y direction
n_no = n_no_x * n_no_y
n_no_active = (n_no_x - 1) * n_no_y #excluding Dirichlet nodes
ele = ti.Vector(4, dt = ti.i32, shape = (n_no_y - 1, n_no_x - 1)) #connectivity matrix
init_pos = ti.var(dt = ti.f32, shape = n_no * 2) # store initial node position (not updated for infinitesimal deformation)
@ti.func
def k_ele(x0, x1, x2, x3):
l_x = ti.Vector([-1, 1, 1, -1])
l_y = ti.Vector([-1, -1, 1, 1])
k = ti.Matrix([[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0]])
a = (x1[0] - x0[0]) / 2
b = (x2[1] - x0[1]) / 2
for i, j in ti.ndrange(4, 4):
# Here is the problem:
k[2 * i, 2 * j] = b / a * (1 + 1 / 3 * l_y[i] * l_y[j]) * l_x[i] * l_x[j] + a / b * (1 - nu) / 2 * (1 + l_x[i] * l_x[j] / 3) * l_y[i] * l_y[j]
return k
@ti.kernel
def init():
for i, j in ele:
ele[i, j][0] = j + i * n_no_x
ele[i, j][1] = ele[i, j][0] + 1
ele[i, j][2] = ele[i, j][1] + n_no_x
ele[i, j][3] = ele[i, j][2] - 1
for i, j in ti.ndrange(n_no_y, n_no_x):
init_pos[2 * (j + i * n_no_x)] = j * (l / (n_no_x - 1))
init_pos[2 * (j + i * n_no_x) + 1] = i * (l / (n_no_x - 1))
for i, j in ele:
x0 = ti.Vector([init_pos[2 * ele[i, j][0]], init_pos[2 * ele[i, j][0] + 1]])
x1 = ti.Vector([init_pos[2 * ele[i, j][1]], init_pos[2 * ele[i, j][1] + 1]])
x2 = ti.Vector([init_pos[2 * ele[i, j][2]], init_pos[2 * ele[i, j][2] + 1]])
x3 = ti.Vector([init_pos[2 * ele[i, j][3]], init_pos[2 * ele[i, j][3] + 1]])
k = k_ele(x0, x1, x2, x3)
init()