在 Matplotlib.pyplot 中创建一个与特征空间垂直的立方体

Question

我正在尝试制作一个所有边都垂直于每个特征向量的立方体，作为一种在 3d 中给定任何可能的法向应力和剪切应力的情况下可视化主应力的方法。我试过使用简单的旋转矩阵并将它们应用于点列表，但似乎总是存在一些错误，我不确定它是否与我应用旋转矩阵的方式、我给它们的角度或顺序有关我用。

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()

我认为最有可能出现错误的地方是我使用的角度、变量 x、y 和 z（以弧度为单位）或我应用旋转矩阵的方法。然而，我可能会完全错误地解决这个问题，也许我应该使用其他一些旋转点的方法。

我的主要目标是将立方体与特征向量对齐，因为立方体只是一种视觉表示。如果您有任何想法或知道我可以做到这一点的不同方法，我们将不胜感激！

Answer 1

我弄明白了，所以基本上角度是错误的，旋转矩阵不够通用。

我首先需要做的是将立方体与其中一个轴对齐，我通过仅围绕 z 轴旋转它到我的一个向量，然后沿着垂直于第一个向量的新向量旋转它来做到这一点xy 平面，最后沿着矢量本身。这是我用来解决这个问题的代码：（我还将添加一个 link 到带有旋转矩阵的 wikapedia 文章，因为那是我想出来的）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)

在 Matplotlib.pyplot 中创建一个与特征空间垂直的立方体

Creating a cube that is normal to an eigenspace in Matplotlib.pyplot

python

matplotlib

matrix

rotation

eigenvector