Rotation matrices don't match reference solution? about continuum_mechanics

[leon-dae]

Update

I dug a little bit deeper.

I compared the loop solution with an updated modulo function that takes into account 0-start-indexing and a written out matrix from p.89 "M. A. Slawinski, Seismic waves and rays in elastic media. Amsterdam: Pergamon, 2003." It does not matter which rotation I make, the looped and fully typed rotation matrix are the same and as such don't match the reference solution in the same places.

I also compared to other reference solutions, e.g. this one.

Something is not correct and I am not able to find the error.

theta = sym.pi/6
c = sym.cos(theta)
s = sym.sin(theta)
Q = sym.Matrix([[1, 0, 0], [0, c, -s], [0, s, c]])      # counterlockwise, x-axis, 30°


def mod3(k):
    k_ = k+1
    if k_ <= 3:
        return k
    return k-3


M = sym.zeros(6)
for i in range(3):
    for j in range(3):
        M[i, j] = Q[i, j]**2
        M[i, j + 3] = 2 * Q[i, mod3(j + 1)] * Q[i, mod3(j + 2)]
        M[i + 3, j] = Q[mod3(i + 1), j] * Q[mod3(i + 2), j]
        M[i + 3, j + 3] = Q[mod3(i + 1), mod3(j + 1)] * \
            Q[mod3(i + 2), mod3(j + 2)] + \
            Q[mod3(i + 1), mod3(j + 2)] * \
            Q[mod3(i + 2), mod3(j + 1)]

M_full = np.array([[Q[0, 0]**2, Q[0, 1]**2, Q[0, 2]**2, 2*Q[0, 1]*Q[0, 2], 2*Q[0, 0]*Q[0, 2], 2*Q[0, 0]*Q[0, 1]],
                   [Q[1, 0]**2, Q[1, 1]**2, Q[1, 2]**2, 2*Q[1, 1]*Q[1, 2], 2*Q[1, 0]*Q[1, 2], 2*Q[1, 0]*Q[1, 1]],
                   [Q[2, 0]**2, Q[2, 1]**2, Q[2, 2]**2, 2*Q[2, 1]*Q[2, 2], 2*Q[2, 0]*Q[2, 2], 2*Q[2, 0]*Q[2, 1]],
                   [Q[1, 0]*Q[2, 0], Q[1, 1]*Q[2, 1], Q[1, 2]*Q[2, 2], Q[1, 1]*Q[2, 2] + Q[1, 2]*Q[2, 1], Q[1, 0]*Q[2, 2] + Q[1, 2]*Q[2, 0], Q[1, 0]*Q[2, 1] + Q[1, 1]*Q[2, 0]],
                   [Q[0, 0]*Q[2, 0], Q[0, 1]*Q[2, 1], Q[0, 2]*Q[2, 2], Q[0, 1]*Q[2, 2] + Q[0, 2]*Q[2, 1], Q[0, 0]*Q[2, 2] + Q[0, 2]*Q[2, 0], Q[0, 0]*Q[2, 1] + Q[0, 1]*Q[2, 1]],
                   [Q[0, 0]*Q[1, 0], Q[0, 1]*Q[1, 1], Q[0, 2]*Q[1, 2], Q[0, 1]*Q[1, 2] + Q[0, 2]*Q[1, 1], Q[0, 0]*Q[1, 2] + Q[0, 2]*Q[1, 0], Q[0, 0]*Q[1, 1] + Q[0, 1]*Q[1, 0]]])

# Counterlockwise, x-axis
M_ref = np.array([[1, 0, 0, 0, 0, 0],
                  [0, c**2, s**2, 2*c*s, 0, 0],
                  [0, s**2, c**2, -2*c*s, 0, 0],
                  [0, -c*s, c*s, c**2-s**2, 0, 0],
                  [0, 0, 0, 0, c, -s],
                  [0, 0, 0, 0, s, c]])

print(M_ref == np.array(M))
print(M_ref == M_full)
print((M_ref == np.array(M)) == (M_ref == M_full))

[leon-dae]

I solved it. The problem was about active and passive rotation. I used the same angle for both the reference and the computed rotation. It appears that for the reference the sign of the angle had to be reversed.

See also the corresponding physics stackexchange post