I am trying to make a cube where all the sides are normal to each of the eigenvectors, as a way to visualize principle stresses given any possible normal and shear stresses in 3d. I've tried using simple rotation matrices and applying them to a list of points but there always seems to be some error and I'm not sure if it is in how I apply the rotation matrices, the angles I give them, or the order I use.
import numpy as np
import math
import matplotlib.pyplot as plt
from math import cos, sin
ax = plt.axes(projection="3d")
# normal and shear stresses
Fx = float(input("Sigma X: "))
Fy = float(input("Sigma Y: "))
Fz = float(input("Sigma Z: "))
Vxy = float(input("Tau xy: "))
Vxz = float(input("Tau xz: "))
Vyz = float(input("Tau yz: "))
A = [[Fx, Vxy, Vxz],
[Vxy, Fy, Vyz],
[Vxz, Vyz, Fz]]
eigval, eigvect = np.linalg.eig(A) # finding principle stresses and their directions
# rounding off error
eigvect = np.round(eigvect, 5)
eigval = np.round (eigval, 5)
# drawing eigenvectors or principle force directions
ax.quiver(0, 0, 0, eigvect[0, 0] * 2, eigvect[1, 0] * 2, eigvect[2, 0] * 2, color="orange")
ax.quiver(0, 0, 0, eigvect[0, 1] * 2, eigvect[1, 1] * 2, eigvect[2, 1] * 2, color="blue")
ax.quiver(0, 0, 0, eigvect[0, 2] * 2, eigvect[1, 2] * 2, eigvect[2, 2] * 2, color="red")
# drawing original normal force directions
ax.quiver(0, 0, 0, 2 * Fx / np.abs(Fx), 0, 0, color="orange", linestyle="dashed")
ax.quiver(0, 0, 0, 0, 2 * Fy / np.abs(Fy), 0, color="blue", linestyle="dashed")
ax.quiver(0, 0, 0, 0, 0, 2 * Fz / np.abs(Fz), color="red", linestyle="dashed")
# points used to draw the cube
points = np.array([[1.0, 1.0, 1.0],
[1.0, -1.0, 1.0],
[-1.0, -1.0, 1.0],
[-1.0, 1.0, 1.0],
[1.0, 1.0, 1.0],
[1.0, 1.0, -1.0],
[1.0, -1.0, -1.0],
[1.0, -1.0, 1.0],
[1.0, -1.0, -1.0],
[-1.0, -1.0, -1.0],
[-1.0, -1.0, 1.0],
[-1.0, -1.0, -1.0],
[-1.0, 1.0, -1.0],
[-1.0, 1.0, 1.0],
[-1.0, 1.0, -1.0],
[1.0, 1.0, -1.0]])
# finding the angles that i need to rotate the cube
x = np.arctan2(eigvect[1, 2], eigvect[2, 2])
y = np.arctan2(eigvect[2, 0], eigvect[0, 0])
z = np.arctan2(eigvect[0, 1], eigvect[1, 1])
# rotation matrices
xrot = np.matrix([[1, 0, 0], [0, cos(y), -sin(y)], [0, sin(y), cos(y)]])
yrot = np.matrix([[cos(x), 0, sin(x)], [0, 1, 0], [-sin(x), 0, cos(x)]])
zrot = np.matrix([[cos(z), -sin(z), 0], [sin(z), cos(z), 0], [0, 0, 1]])
# updating the points list to rotate the cube
for i in range(len(points)):
points[i] = np.dot(points[i], xrot)
points[i] = np.dot(points[i], yrot)
points[i] = np.dot(points[i], zrot)
# plotting the rotated cube
ax.plot(points[:, 0], points[:, 1], points[:, 2], color="blue")
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
plt.show()
The most likely spot for the error is I think in either the angles that im using, variables x, y, and z, (in radians) or in the method I apply the rotation matrices. I could however be going about this completely wrong and maybe I should be using some other method of rotating the points.
My main goal is to just align the cube with the eigenvectors as the cube is just a visual representation. If you have any ideas or know of a different way i could do this it would be greatly appreciated!
I figured it out, so basicaly the angles were wrong and the rotation matrices were not general enough.
What I first needed to do was align the cube with one of the axies, I did this by rotating it only about the z axis to one of my vectors, then rotate it along a new vector perpendicular to the first on the xy plane, and finaly along the vector its self. Heres the code I used to fix this: (I'll also add a link to the wikapedia article with the rotation matrices as that was where i figured it out) https://en.wikipedia.org/wiki/Rotation_matrix
# angles used
theta = - np.arctan2(eigvect[1, 0], eigvect[0, 0])
Vxy = np.arccos(eigvect[2, 0])
Uxyz = - np.arctan2(eigvect[1, 2], eigvect[2, 2])
# z-axis rotation
z_rot = np.matrix([[cos(theta), -sin(theta), 0],
[sin(theta), cos(theta), 0],
[0, 0, 1]])
# making the perpendicular unit vector
magnetude = np.sqrt(eigvect[1, 0] ** 2 + eigvect[0, 0] ** 2)
vx = eigvect[1, 0] / magnetude
vy = - eigvect[0, 0] / magnetude
vz = 0
# xy-plane rotation
Vxy_rot = np.matrix([[cos(Vxy) + (vx ** 2) * (1 - cos(Vxy)), vx * vy * (1 - cos(Vxy)) - vz * sin(Vxy), vx * vz * (1 - cos(Vxy)) + vy * sin(Vxy)],
[vy * vx * (1 - cos(Vxy)) + vz * sin(Vxy), cos(Vxy) + (vy ** 2) * (1 - cos(Vxy)), vy * vz * (1 - cos(Vxy)) - vx * sin(Vxy)],
[vz * vx * (1 - cos(Vxy)) - vy * sin(Vxy), vz * vy * (1 - cos(Vxy)) + vx * sin(Vxy), cos(Vxy) + (vz ** 2) * (1 - cos(Vxy))]])
# x, y, z components
ux = eigvect[0, 0]
uy = eigvect[1, 0]
uz = eigvect[2, 0]
# xyz rotation
Uxyz_rot = np.matrix([[cos(Uxyz) + (ux ** 2) * (1 - cos(Uxyz)), ux * uy * (1 - cos(Uxyz)) - uz * sin(Uxyz), ux * uz * (1 - cos(Uxyz)) + uy * sin(Uxyz)],
[uy * ux * (1 - cos(Uxyz)) + uz * sin(Uxyz), cos(Uxyz) + (uy ** 2) * (1 - cos(Uxyz)), uy * uz * (1 - cos(Uxyz)) - ux * sin(Uxyz)],
[uz * ux * (1 - cos(Uxyz)) - uy * sin(Uxyz), uz * uy * (1 - cos(Uxyz)) + ux * sin(Uxyz), cos(Uxyz) + (uz ** 2) * (1 - cos(Uxyz))]])
# applying the rotations
for i in range(len(points)):
points[i] = np.dot(points[i], z_rot)
points[i] = np.dot(points[i], Vxy_rot)
points[i] = np.dot(points[i], Uxyz_rot)