Skip to content
Snippets Groups Projects
Commit 1183befd authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files

Rotation2

parent d2edcc81
No related branches found
No related tags found
1 merge request!140Rotation2
Pipeline #11742 passed
Stage: test
Stage: deploy
Hide whitespace changes
Inline Side-by-side
wetb/hawc2/st_file.py
+ 109
− 2
View file @ 1183befd
...@@ -38,8 +38,8 @@ class StFile(object): ...@@ -38,8 +38,8 @@ class StFile(object):
- Ix : area moment of inertia with respect to principal bending xe axis [m4]. This is the principal bending axis most parallel to the xc2 axis - Ix : area moment of inertia with respect to principal bending xe axis [m4]. This is the principal bending axis most parallel to the xc2 axis
- Iy : area moment of inertia with respect to principal bending ye axis [m4] - Iy : area moment of inertia with respect to principal bending ye axis [m4]
- K : torsional stiffness constant with respect to ze axis at the shear center [m4/rad]. For a circular section only this is identical to the polar moment of inertia. - K : torsional stiffness constant with respect to ze axis at the shear center [m4/rad]. For a circular section only this is identical to the polar moment of inertia.
- kx : shear factor for force in principal bending xe direction [-] - k_x : shear factor for force in principal bending xe direction [-]
- ky : shear factor for force in principal bending ye direction [-] - k_y : shear factor for force in principal bending ye direction [-]
- A : cross sectional area [m2] - A : cross sectional area [m2]
- pitch : structural pitch about z_c2 axis. This is the angle between the xc2 -axis defined with the c2_def command and the main principal bending axis xe. - pitch : structural pitch about z_c2 axis. This is the angle between the xc2 -axis defined with the c2_def command and the main principal bending axis xe.
- xe : xc2-coordinate from C1/2 to center of elasticity [m]. The elastic center is the point where radial force (in the z-direction) does not contribute to bending around the x or y directions. - xe : xc2-coordinate from C1/2 to center of elasticity [m]. The elastic center is the point where radial force (in the z-direction) does not contribute to bending around the x or y directions.
...@@ -142,6 +142,113 @@ class StFile(object): ...@@ -142,6 +142,113 @@ class StFile(object):
fid.write('$%i %i\n' % (set, npoints)) fid.write('$%i %i\n' % (set, npoints))
fid.write(dstr + '\n') fid.write(dstr + '\n')
def element_stiffnessmatrix(self, radius, mset_nr, set_nr, length):
"""Compute the element stiffness matrix
Parameters
----------
radius : float
radius of element (used of obtain element properties
length : float
eleement length
"""
K = np.zeros((13, 13))
"r m x_cg y_cg ri_x ri_y x_sh y_sh E G I_x I_y I_p k_x k_y A pitch x_e y_e"
ES1, ES2, EMOD, GMOD, IX, IY, IZ, KX, KY, A = [getattr(self, n)(radius, mset_nr, set_nr)
for n in "x_sh,y_sh,E,G,I_x,I_y,I_p,k_x,k_y,A".split(",")]
ELLGTH = length
ETAX = EMOD * IX / (KY * GMOD * A * ELLGTH**2)
ETAY = EMOD * IY / (KX * GMOD * A * ELLGTH**2)
ROX = 1 / (1 + 12 * ETAX)
ROY = 1 / (1 + 12 * ETAY)
ROY = 1 / (1 + 12 * ETAX)
K[1, 1] = 12 * EMOD * IY * ROY / ELLGTH**3
K[1, 5] = 6 * EMOD * IY * ROY / ELLGTH**2
K[1, 6] = -K[1, 1] * ES2
K[1, 7] = -K[1, 1]
K[1, 11] = K[1, 5]
K[1, 12] = -K[1, 6]
K[2, 2] = 12 * EMOD * IX * ROX / ELLGTH**3
K[2, 4] = -6 * EMOD * IX * ROX / ELLGTH**2
K[2, 6] = K[2, 2] * ES1
K[2, 8] = -K[2, 2]
K[2, 10] = K[2, 4]
K[2, 12] = -K[2, 6]
K[3, 3] = A * EMOD / ELLGTH
K[3, 9] = -K[3, 3]
K[4, 4] = 4 * EMOD * IX * (1 + 3 * ETAX) * ROX / ELLGTH
K[4, 6] = K[2, 4] * ES1
K[4, 8] = -K[2, 4]
K[4, 10] = 2 * EMOD * IX * (1 - 6 * ETAX) * ROX / ELLGTH
K[4, 12] = -K[4, 6]
K[5, 5] = 4 * EMOD * IY * (1 + 3 * ETAY) * ROY / ELLGTH
K[5, 6] = -K[1, 5] * ES2
K[5, 7] = -K[1, 5]
K[5, 11] = 2 * EMOD * IY * (1 - 6 * ETAY) * ROY / ELLGTH
K[5, 12] = -K[5, 6]
K[6, 6] = GMOD * IZ / ELLGTH + 12 * EMOD * (IX * ES1**2 * ROX + IY * ES2**2 * ROY) / ELLGTH**3
K[6, 7] = K[1, 12]
K[6, 8] = K[2, 12]
K[6, 10] = -K[4, 12]
K[6, 11] = -K[5, 12]
K[6, 12] = -K[6, 6]
K[7, 7] = K[1, 1]
K[7, 11] = -K[1, 5]
K[7, 12] = K[1, 6]
K[8, 8] = K[2, 2]
K[8, 10] = -K[2, 4]
K[8, 12] = K[2, 6]
K[9, 9] = K[3, 3]
K[10, 10] = K[4, 4]
K[10, 12] = -K[4, 6]
K[11, 11] = K[5, 5]
K[11, 12] = -K[5, 6]
K[12, 12] = K[6, 6]
K = K[1:, 1:]
K = K + K.T - np.eye(12) * K
return K
def shape_function_ori(self, radius, mset_nr, set_nr, length, z):
XSC, YSC, EMOD, GMOD, IX, IY, IZ, KX, KY, AREA = [getattr(self, n)(radius, mset_nr, set_nr)
for n in "x_sh,y_sh,E,G,I_x,I_y,I_p,k_x,k_y,A".split(",")]
etax = EMOD * IX / KY / GMOD / AREA / (length**2)
etay = EMOD * IY / KX / GMOD / AREA / (length**2)
rhox = 1 / (1 + 12 * etax)
rhoy = 1 / (1 + 12 * etay)
f1 = z / length
f2 = 1 - f1
f3x = f1 * (3 * f1 - 2 * f1**2 + 12 * etax) * rhox
f4x = f2 * (3 * f2 - 2 * f2**2 + 12 * etax) * rhox
f5x = length * (f1**2 - (1 - 6 * etax) * f1 - 6 * etax) * f1 * rhox
f6x = length * (f2**2 - (1 - 6 * etax) * f2 - 6 * etax) * f2 * rhox
f7x = 6 / length * f1 * f2 * rhox
f8x = f1 * (3 * f1 - 2 * (1 - 6 * etax)) * rhox
f9x = f2 * (3 * f2 - 2 * (1 - 6 * etax)) * rhox
f3y = f1 * (3 * f1 - 2 * f1**2 + 12 * etay) * rhoy
f4y = f2 * (3 * f2 - 2 * f2**2 + 12 * etay) * rhoy
f5y = length * (f1**2 - (1 - 6 * etay) * f1 - 6 * etay) * f1 * rhoy
f6y = length * (f2**2 - (1 - 6 * etay) * f2 - 6 * etay) * f2 * rhoy
f7y = 6 / length * f1 * f2 * rhoy
f8y = f1 * (3 * f1 - 2 * (1 - 6 * etay)) * rhoy
f9y = f2 * (3 * f2 - 2 * (1 - 6 * etay)) * rhoy
return np.array([[f4y, 0, 0, 0, -f7y, 0],
[0, f4x, 0, f7x, 0, 0],
[0, 0, f2, 0, 0, 0],
[0, f6x, 0, f9x, 0, 0],
[-f6y, 0, 0, 0, f9y, 0],
[(f2 - f4y) * YSC, -(f2 - f4x) * XSC, 0, f7x * XSC, f7y * YSC, f2],
[f3y, 0, 0, 0, f7y, 0],
[0, f3x, 0, -f7x, 0, 0],
[0, 0, f1, 0, 0, 0],
[0, -f5x, 0, f8x, 0, 0],
[f5y, 0, 0, 0, f8y, 0],
[(f1 - f3y) * YSC, -(f1 - f3x) * XSC, 0, -f7x * XSC, -f7y * YSC, f1]]).T
if __name__ == "__main__": if __name__ == "__main__":
import os import os
......
wetb/utils/rotation.py
+ 52
− 8
View file @ 1183befd
...@@ -6,6 +6,28 @@ Created on 21/07/2014 ...@@ -6,6 +6,28 @@ Created on 21/07/2014
import numpy as np import numpy as np
def s2matrix(s):
"""Creates a transformation matrix from string
Parameters
----------
s : string
x is
Returns
-------
matrix : array_like
3x3 transformation matrix
Examples
--------
>> s2matrix('x,y,z') # identity matrix
>> s2matrix('x,-z,y) # 90 deg rotation around x
"""
d = {xyz: v for xyz, v in zip('xyz', np.eye(3))}
return np.array([eval(xyz, d) for xyz in s.split(",")]).T
def transformation_matrix(angles, xyz): def transformation_matrix(angles, xyz):
"""Create Transformation matrix(es) """Create Transformation matrix(es)
!!!Note that the columns of the returned matrix(es) is original(unrotate) xyz-axes in rotated coordinates\n !!!Note that the columns of the returned matrix(es) is original(unrotate) xyz-axes in rotated coordinates\n
...@@ -266,24 +288,34 @@ def norm(vector): ...@@ -266,24 +288,34 @@ def norm(vector):
# axis to ... # axis to ...
#======================================================================================================================= #=======================================================================================================================
def axis2axis_angle(axis): def axis2axis_angle(axis):
"""
deg : boolean
if True, axis length is assumed to be angle in deg
"""
axis = np.asarray(axis) axis = np.asarray(axis)
angle = np.sqrt(((axis**2).sum())) angle = np.sqrt(((axis**2).sum()))
x, y, z = axis / np.sqrt((axis**2).sum()) return np.r_[axis / np.sqrt((axis**2).sum()), angle]
return np.r_[axis / np.sqrt((axis**2).sum()), np.deg2rad(angle)]
def axis2matrix(axis): def axis2matrix(axis, deg=False):
# http://www.euclideanspace.com/maths/geometry/rotations/conversions/angleToMatrix/index.htm
axis = np.asarray(axis) axis = np.asarray(axis)
angle = np.sqrt(((axis**2).sum())) angle = np.sqrt(((axis**2).sum()))
x, y, z = axis / np.sqrt((axis**2).sum()) if deg:
angle = np.deg2rad(angle)
if angle == 0:
return np.eye(3)
x, y, z = xyz = axis / np.sqrt((axis**2).sum())
angle = np.deg2rad(angle)
c, s = np.cos(angle), np.sin(angle) c, s = np.cos(angle), np.sin(angle)
t = 1 - c t = 1 - c
return np.array([[t * x * x + c, t * x * y - z * s, t * x * z + y * s], return np.array([[t * x * x + c, t * x * y - z * s, t * x * z + y * s],
[t * x * y + z * s, t * y * y + c, t * y * z - x * s], [t * x * y + z * s, t * y * y + c, t * y * z - x * s],
[t * x * z - y * s, t * y * z + x * s, t * z * z + c]]) [t * x * z - y * s, t * y * z + x * s, t * z * z + c]])
# alternative implementation
#asym = np.array([[0, -z, y], [z, 0, -x], [-y, x, 0]])
# return c * np.eye(3, 3) + s * asym + (1 - c) * xyz[np.newaxis] * xyz[:, np.newaxis]
#======================================================================================================================= #=======================================================================================================================
...@@ -293,11 +325,13 @@ def axis2matrix(axis): ...@@ -293,11 +325,13 @@ def axis2matrix(axis):
def axis_angle2axis(axis_angle): def axis_angle2axis(axis_angle):
x, y, z, angle = axis_angle x, y, z, angle = axis_angle
return np.array([x, y, z]) * np.rad2deg(angle) return np.array([x, y, z]) * angle
def axis_angle2quaternion(axis_angle): def axis_angle2quaternion(axis_angle, deg=False):
x, y, z, angle = axis_angle x, y, z, angle = axis_angle
if deg:
angle = np.deg2rad(angle)
s = np.sin(angle / 2) s = np.sin(angle / 2)
x = x * s x = x * s
y = y * s y = y * s
...@@ -310,9 +344,11 @@ def axis_angle2quaternion(axis_angle): ...@@ -310,9 +344,11 @@ def axis_angle2quaternion(axis_angle):
# # quaternion to ... # # quaternion to ...
#======================================================================================================================= #=======================================================================================================================
def quaternion2axis_angle(quaternion): def quaternion2axis_angle(quaternion, deg=False):
qw, qx, qy, qz = quaternion / norm(quaternion) qw, qx, qy, qz = quaternion / norm(quaternion)
angle = 2 * np.arccos(qw) angle = 2 * np.arccos(qw)
if deg:
angle = np.rad2deg(angle)
t = np.sqrt(1 - qw**2) t = np.sqrt(1 - qw**2)
x = qx / t x = qx / t
y = qy / t y = qy / t
...@@ -373,3 +409,11 @@ def matrix2quaternion(matrix): ...@@ -373,3 +409,11 @@ def matrix2quaternion(matrix):
qz = 0.25 * S qz = 0.25 * S
e = np.array([qw, qx, qy, qz]) e = np.array([qw, qx, qy, qz])
return e / np.sqrt(np.sum(e**2)) return e / np.sqrt(np.sum(e**2))
def matrix2axis_angle(matrix, deg=False):
return quaternion2axis_angle(matrix2quaternion(matrix), deg)
def matrix2axis(matrix, deg=False):
return axis_angle2axis(matrix2axis_angle(matrix, deg))
wetb/utils/tests/test_rotation.py
+ 9
− 4
View file @ 1183befd
...@@ -9,13 +9,18 @@ import numpy as np ...@@ -9,13 +9,18 @@ import numpy as np
from wetb.utils.geometry import rad from wetb.utils.geometry import rad
from wetb.utils.rotation import transformation_matrix, mdot, dots, rotate, rotate_x, \ from wetb.utils.rotation import transformation_matrix, mdot, dots, rotate, rotate_x, \
rotmat, rotate_y, rotate_z, norm, axis2axis_angle, axis_angle2axis, axis2matrix, axis_angle2quaternion,\ rotmat, rotate_y, rotate_z, norm, axis2axis_angle, axis_angle2axis, axis2matrix, axis_angle2quaternion,\
quaternion2matrix, quaternion2axis_angle, matrix2quaternion quaternion2matrix, quaternion2axis_angle, matrix2quaternion, s2matrix
from tests import npt from tests import npt
import pytest import pytest
x, y, z = 0, 1, 2 x, y, z = 0, 1, 2
def test_s2matrix():
npt.assert_array_equal(s2matrix('x,-y,-z'), np.array([[1, 0, 0], [0, -1, 0], [0, 0, -1]]).T)
npt.assert_array_equal(s2matrix('x,-z,y'), np.array([[1, 0, 0], [0, 0, -1], [0, 1, 0]]).T)
def test_transformation_matrix(): def test_transformation_matrix():
npt.assert_array_almost_equal(transformation_matrix(rad(30), x), [[[1., 0., 0.], npt.assert_array_almost_equal(transformation_matrix(rad(30), x), [[[1., 0., 0.],
[0., 0.8660254, 0.5], [0., 0.8660254, 0.5],
...@@ -136,7 +141,7 @@ def test_norm(): ...@@ -136,7 +141,7 @@ def test_norm():
def test_axis2axis_angle(): def test_axis2axis_angle():
axis_angle = np.array([2.504687681842697E-002, 0.654193239950699, 0.755912599950848, 3.09150937138433]) axis_angle = np.array([2.504687681842697E-002, 0.654193239950699, 0.755912599950848, np.rad2deg(3.09150937138433)])
axis = np.array([4.43656429408, 115.877535983, 133.895130906]) axis = np.array([4.43656429408, 115.877535983, 133.895130906])
npt.assert_array_almost_equal(axis2axis_angle(axis), axis_angle) npt.assert_array_almost_equal(axis2axis_angle(axis), axis_angle)
npt.assert_array_almost_equal(axis_angle2axis(axis2axis_angle(axis)), axis) npt.assert_array_almost_equal(axis_angle2axis(axis2axis_angle(axis)), axis)
...@@ -156,11 +161,11 @@ axis_quaternion_matrix_lst = [([30, 0, 0], [0.9659258262890683, 0.25881904510252 ...@@ -156,11 +161,11 @@ axis_quaternion_matrix_lst = [([30, 0, 0], [0.9659258262890683, 0.25881904510252
@pytest.mark.parametrize("axis,quaternion,matrix", axis_quaternion_matrix_lst) @pytest.mark.parametrize("axis,quaternion,matrix", axis_quaternion_matrix_lst)
def test_axis2matrix(axis, quaternion, matrix): def test_axis2matrix(axis, quaternion, matrix):
npt.assert_array_almost_equal(axis2matrix(axis), matrix) npt.assert_array_almost_equal(axis2matrix(axis, deg=True), matrix)
def test_axis_angle2axis(): def test_axis_angle2axis():
axis_angle = np.array([2.504687681842697E-002, 0.654193239950699, 0.755912599950848, 3.09150937138433]) axis_angle = np.array([2.504687681842697E-002, 0.654193239950699, 0.755912599950848, np.rad2deg(3.09150937138433)])
axis = np.array([4.43656429408, 115.877535983, 133.895130906]) axis = np.array([4.43656429408, 115.877535983, 133.895130906])
npt.assert_array_almost_equal(axis_angle2axis(axis_angle), axis) npt.assert_array_almost_equal(axis_angle2axis(axis_angle), axis)
npt.assert_array_almost_equal(axis2axis_angle(axis_angle2axis(axis_angle)), axis_angle) npt.assert_array_almost_equal(axis2axis_angle(axis_angle2axis(axis_angle)), axis_angle)
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment