关于ti.sym_eig

综合讨论区
#1 
上下cauchy_stress = ti.Matrix([[0.014860193, 0.0], [0.0, 0.014860193]])值相同

直接给出

    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)

cauchy_stress矩阵是相同的值，想进行特征分解，给定矩阵值能正确分解。

不直接赋值矩阵，就分解错误。也没有报错，请问这是什么原因。

不直接给出，

    # 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)
#2

不好意思没太理解你的问题，如果不直接给出cauchy_stress的话，if all(cauchy_stress.transpose() == cauchy_stress):这一步要怎么比较呢 :thinking:

#3

本身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就可以分解正确，得到特征向量是单位向量。所以就不太明白为什么会这样。

#4

Hi @wjy. 欢迎来到太极论坛。

不太清楚你的问题。
你的问题是通过计算出来的矩阵直接做特征分解和常量矩阵的特征分解结果不同？

#5

对的，但是常量矩阵和计算出来的矩阵是相等的

#8

源码放在这里了，


实在不知道哪里的问题了，请您帮忙康康，谢谢 :sob:
#9

抱歉，我没有石墨文档的账号。能贴个gist或者github链接么？

PS：最好能简化一下代码，到最小的复现你问题的程度。目前来看代码有点多，开发者看起来会很费劲。

#10

https://github.com/JiaYanWang-Home/damage
放在了这里，简化了，您康康 :grinning: