Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Rotation matrices don't match reference solution? #178

Closed
leon-dae opened this issue Oct 24, 2020 · 2 comments
Closed

Rotation matrices don't match reference solution? #178

leon-dae opened this issue Oct 24, 2020 · 2 comments

Comments

@leon-dae
Copy link

  • Continuum Mechanics version: 0.3.0

Description

The rotation matrix given in http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php under section 3.2.11 is not the same as what I get from the package. It might be an issue with clockwise and counterclockwise rotation, as the numbers are the same, the signs are not always.

As far as I understand it, a counterclockwise rotation should be positive about the rotation axis according to the right-hand rule. This is also what's shown in his book and the figure under 3.2.11 on the webpage.

I know this is a version under development, but I checked your code and did not find any error.

What I Did

Computed the two rotation matrices as below. They don't match.

theta = sym.pi/6
c = sym.cos(theta)
s = sym.sin(theta)

# counterclockwise, z-axis, 30°
# Taken from https://en.wikipedia.org/wiki/Rotation_matrix#Basic_rotations
Q = sym.Matrix([[c, -s, 0], [s, c, 0], [0, 0, 1]])     

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, (j + 1) % 3] * Q[i, (j + 2) % 3]
        M[i + 3, j] = Q[(i + 1) % 3, j] * Q[(i + 2) % 3, j]
        M[i + 3, j + 3] = Q[(i + 1) % 3, (j + 1) % 3] * \
            Q[(i + 2) % 3, (j + 2) % 3] + \
            Q[(i + 1) % 3, (j + 2) % 3] * \
            Q[(i + 2) % 3, (j + 1) % 3]

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

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

Output:

[[ True  True  True  True  True False]
 [ True  True  True  True  True False]
 [ True  True  True  True  True  True]
 [ True  True  True  True  True  True]
 [ True  True  True  True  True  True]
 [False False  True  True  True  True]]
@leon-dae
Copy link
Author

leon-dae commented Oct 27, 2020

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
Copy link
Author

leon-dae commented Nov 5, 2020

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

@leon-dae leon-dae closed this as completed Nov 5, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant