本身cauchy_stress是通过每一步计算得到的
@ti.kernel
def update_d():
# d映射到网格 d:P2G
map_d_to_grid()
# 计算 d 的拉普拉斯算子
central_laplacian()
for p in range(particle_num[None]):
xp = particles[p].position / dx
base = int(xp - 0.5)
fx = xp - base
w = [0.5 * (1.5 - fx) ** 2, 0.75 - (fx - 1) ** 2, 0.5 * (fx - 0.5) ** 2]
# 在网格上计算d的拉普拉斯算子之后,映射回粒子,得到粒子上的拉普拉斯算子 p_laplacian_d
p_laplacian_d = 0.0
for i, j in ti.static(ti.ndrange(3, 3)): # loop over 3x3 grid node neighborhood
g_laplacian_d = node[base + ti.Vector([i, j])].laplacian_d
weight = w[i][0] * w[j][1]
p_laplacian_d += weight * g_laplacian_d
# 计算阻抗 compute geometric resistance:Dc
Dc = particles[p].d - l0 ** 2 * p_laplacian_d
# print("Dc:", Dc)
# *********计算未退化(g(d) = 1)的柯西应力 compute un-dearaded Cauchy stress***********
# ++++++++ QR分解和符号规定 +++++++++++
Q, R = QR2(particles[p].F)
# sign convention 符号规定
# 之前直接Q[0, 0]报错,取不到值
r00 = R[0, 0]
if r00 < 0:
Q *= -1
R *= -1
# 计算 J
particles[p].J = 1.0
for r in ti.static(range(2)):
particles[p].J *= R[r, r]
# 计算 stress
kx = 0 * mu_0
pk1 = calculate_pk1(Q, R, particles[p].J, kx, 1.0)
# P = Jσ(F-T)(F的逆的转置) σ = PF/J
cauchy_stress = ti.Matrix.zero(float, 2, 2)
cauchy_stress = pk1 @ particles[p].F / particles[p].J
# print(cauchy_stress)
# ******************************
# 计算sigma plus
sigma_plus = ti.Matrix.zero(float, 2, 2)
T1 = ti.Matrix.zero(float, 2, 2)
T2 = ti.Matrix.zero(float, 4, 2)
sym_eigT1 = ti.Vector.zero(float, 2)
sym_eigT2 = ti.Matrix.zero(float, 2, 2)
# 直接给出
# cauchy_stress = ti.Matrix([[0.014860193, 0.0], [0.0, 0.014860193]])
# 计算特征值和特征矩阵
if all(cauchy_stress.transpose() == cauchy_stress):
sym_eigT1, sym_eigT2 = ti.sym_eig(cauchy_stress, ti.f32)
eigen_values = sym_eigT1
eigen_vector = sym_eigT2
print("cauchy_stress:", cauchy_stress, "eig1:", sym_eigT1, "eig2:", sym_eigT2)
else:
T1,T2 = ti.eig(cauchy_stress, ti.f32)
eigen_values = T1
eigen_vector = T2
其中一步得到的值输出看了一下是[[0.014860193, 0.0], [0.0, 0.014860193]],用sym_eig分解得到的不是单位向量,测试的时候直接给cauchy_stress赋值[[0.014860193, 0.0], [0.0, 0.014860193]]同样用sym_eig就可以分解正确,得到特征向量是单位向量。所以就不太明白为什么会这样。