test_body.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340
  1. from sympy import (Symbol, symbols, sin, cos, Matrix, zeros,
  2. simplify)
  3. from sympy.physics.vector import Point, ReferenceFrame, dynamicsymbols, Dyadic
  4. from sympy.physics.mechanics import inertia, Body
  5. from sympy.testing.pytest import raises, warns_deprecated_sympy
  6. def test_default():
  7. with warns_deprecated_sympy():
  8. body = Body('body')
  9. assert body.name == 'body'
  10. assert body.loads == []
  11. point = Point('body_masscenter')
  12. point.set_vel(body.frame, 0)
  13. com = body.masscenter
  14. frame = body.frame
  15. assert com.vel(frame) == point.vel(frame)
  16. assert body.mass == Symbol('body_mass')
  17. ixx, iyy, izz = symbols('body_ixx body_iyy body_izz')
  18. ixy, iyz, izx = symbols('body_ixy body_iyz body_izx')
  19. assert body.inertia == (inertia(body.frame, ixx, iyy, izz, ixy, iyz, izx),
  20. body.masscenter)
  21. def test_custom_rigid_body():
  22. # Body with RigidBody.
  23. rigidbody_masscenter = Point('rigidbody_masscenter')
  24. rigidbody_mass = Symbol('rigidbody_mass')
  25. rigidbody_frame = ReferenceFrame('rigidbody_frame')
  26. body_inertia = inertia(rigidbody_frame, 1, 0, 0)
  27. with warns_deprecated_sympy():
  28. rigid_body = Body('rigidbody_body', rigidbody_masscenter,
  29. rigidbody_mass, rigidbody_frame, body_inertia)
  30. com = rigid_body.masscenter
  31. frame = rigid_body.frame
  32. rigidbody_masscenter.set_vel(rigidbody_frame, 0)
  33. assert com.vel(frame) == rigidbody_masscenter.vel(frame)
  34. assert com.pos_from(com) == rigidbody_masscenter.pos_from(com)
  35. assert rigid_body.mass == rigidbody_mass
  36. assert rigid_body.inertia == (body_inertia, rigidbody_masscenter)
  37. assert rigid_body.is_rigidbody
  38. assert hasattr(rigid_body, 'masscenter')
  39. assert hasattr(rigid_body, 'mass')
  40. assert hasattr(rigid_body, 'frame')
  41. assert hasattr(rigid_body, 'inertia')
  42. def test_particle_body():
  43. # Body with Particle
  44. particle_masscenter = Point('particle_masscenter')
  45. particle_mass = Symbol('particle_mass')
  46. particle_frame = ReferenceFrame('particle_frame')
  47. with warns_deprecated_sympy():
  48. particle_body = Body('particle_body', particle_masscenter,
  49. particle_mass, particle_frame)
  50. com = particle_body.masscenter
  51. frame = particle_body.frame
  52. particle_masscenter.set_vel(particle_frame, 0)
  53. assert com.vel(frame) == particle_masscenter.vel(frame)
  54. assert com.pos_from(com) == particle_masscenter.pos_from(com)
  55. assert particle_body.mass == particle_mass
  56. assert not hasattr(particle_body, "_inertia")
  57. assert hasattr(particle_body, 'frame')
  58. assert hasattr(particle_body, 'masscenter')
  59. assert hasattr(particle_body, 'mass')
  60. assert particle_body.inertia == (Dyadic(0), particle_body.masscenter)
  61. assert particle_body.central_inertia == Dyadic(0)
  62. assert not particle_body.is_rigidbody
  63. particle_body.central_inertia = inertia(particle_frame, 1, 1, 1)
  64. assert particle_body.central_inertia == inertia(particle_frame, 1, 1, 1)
  65. assert particle_body.is_rigidbody
  66. with warns_deprecated_sympy():
  67. particle_body = Body('particle_body', mass=particle_mass)
  68. assert not particle_body.is_rigidbody
  69. point = particle_body.masscenter.locatenew('point', particle_body.x)
  70. point_inertia = particle_mass * inertia(particle_body.frame, 0, 1, 1)
  71. particle_body.inertia = (point_inertia, point)
  72. assert particle_body.inertia == (point_inertia, point)
  73. assert particle_body.central_inertia == Dyadic(0)
  74. assert particle_body.is_rigidbody
  75. def test_particle_body_add_force():
  76. # Body with Particle
  77. particle_masscenter = Point('particle_masscenter')
  78. particle_mass = Symbol('particle_mass')
  79. particle_frame = ReferenceFrame('particle_frame')
  80. with warns_deprecated_sympy():
  81. particle_body = Body('particle_body', particle_masscenter,
  82. particle_mass, particle_frame)
  83. a = Symbol('a')
  84. force_vector = a * particle_body.frame.x
  85. particle_body.apply_force(force_vector, particle_body.masscenter)
  86. assert len(particle_body.loads) == 1
  87. point = particle_body.masscenter.locatenew(
  88. particle_body._name + '_point0', 0)
  89. point.set_vel(particle_body.frame, 0)
  90. force_point = particle_body.loads[0][0]
  91. frame = particle_body.frame
  92. assert force_point.vel(frame) == point.vel(frame)
  93. assert force_point.pos_from(force_point) == point.pos_from(force_point)
  94. assert particle_body.loads[0][1] == force_vector
  95. def test_body_add_force():
  96. # Body with RigidBody.
  97. rigidbody_masscenter = Point('rigidbody_masscenter')
  98. rigidbody_mass = Symbol('rigidbody_mass')
  99. rigidbody_frame = ReferenceFrame('rigidbody_frame')
  100. body_inertia = inertia(rigidbody_frame, 1, 0, 0)
  101. with warns_deprecated_sympy():
  102. rigid_body = Body('rigidbody_body', rigidbody_masscenter,
  103. rigidbody_mass, rigidbody_frame, body_inertia)
  104. l = Symbol('l')
  105. Fa = Symbol('Fa')
  106. point = rigid_body.masscenter.locatenew(
  107. 'rigidbody_body_point0',
  108. l * rigid_body.frame.x)
  109. point.set_vel(rigid_body.frame, 0)
  110. force_vector = Fa * rigid_body.frame.z
  111. # apply_force with point
  112. rigid_body.apply_force(force_vector, point)
  113. assert len(rigid_body.loads) == 1
  114. force_point = rigid_body.loads[0][0]
  115. frame = rigid_body.frame
  116. assert force_point.vel(frame) == point.vel(frame)
  117. assert force_point.pos_from(force_point) == point.pos_from(force_point)
  118. assert rigid_body.loads[0][1] == force_vector
  119. # apply_force without point
  120. rigid_body.apply_force(force_vector)
  121. assert len(rigid_body.loads) == 2
  122. assert rigid_body.loads[1][1] == force_vector
  123. # passing something else than point
  124. raises(TypeError, lambda: rigid_body.apply_force(force_vector, 0))
  125. raises(TypeError, lambda: rigid_body.apply_force(0))
  126. def test_body_add_torque():
  127. with warns_deprecated_sympy():
  128. body = Body('body')
  129. torque_vector = body.frame.x
  130. body.apply_torque(torque_vector)
  131. assert len(body.loads) == 1
  132. assert body.loads[0] == (body.frame, torque_vector)
  133. raises(TypeError, lambda: body.apply_torque(0))
  134. def test_body_masscenter_vel():
  135. with warns_deprecated_sympy():
  136. A = Body('A')
  137. N = ReferenceFrame('N')
  138. with warns_deprecated_sympy():
  139. B = Body('B', frame=N)
  140. A.masscenter.set_vel(N, N.z)
  141. assert A.masscenter_vel(B) == N.z
  142. assert A.masscenter_vel(N) == N.z
  143. def test_body_ang_vel():
  144. with warns_deprecated_sympy():
  145. A = Body('A')
  146. N = ReferenceFrame('N')
  147. with warns_deprecated_sympy():
  148. B = Body('B', frame=N)
  149. A.frame.set_ang_vel(N, N.y)
  150. assert A.ang_vel_in(B) == N.y
  151. assert B.ang_vel_in(A) == -N.y
  152. assert A.ang_vel_in(N) == N.y
  153. def test_body_dcm():
  154. with warns_deprecated_sympy():
  155. A = Body('A')
  156. B = Body('B')
  157. A.frame.orient_axis(B.frame, B.frame.z, 10)
  158. assert A.dcm(B) == Matrix([[cos(10), sin(10), 0], [-sin(10), cos(10), 0], [0, 0, 1]])
  159. assert A.dcm(B.frame) == Matrix([[cos(10), sin(10), 0], [-sin(10), cos(10), 0], [0, 0, 1]])
  160. def test_body_axis():
  161. N = ReferenceFrame('N')
  162. with warns_deprecated_sympy():
  163. B = Body('B', frame=N)
  164. assert B.x == N.x
  165. assert B.y == N.y
  166. assert B.z == N.z
  167. def test_apply_force_multiple_one_point():
  168. a, b = symbols('a b')
  169. P = Point('P')
  170. with warns_deprecated_sympy():
  171. B = Body('B')
  172. f1 = a*B.x
  173. f2 = b*B.y
  174. B.apply_force(f1, P)
  175. assert B.loads == [(P, f1)]
  176. B.apply_force(f2, P)
  177. assert B.loads == [(P, f1+f2)]
  178. def test_apply_force():
  179. f, g = symbols('f g')
  180. q, x, v1, v2 = dynamicsymbols('q x v1 v2')
  181. P1 = Point('P1')
  182. P2 = Point('P2')
  183. with warns_deprecated_sympy():
  184. B1 = Body('B1')
  185. B2 = Body('B2')
  186. N = ReferenceFrame('N')
  187. P1.set_vel(B1.frame, v1*B1.x)
  188. P2.set_vel(B2.frame, v2*B2.x)
  189. force = f*q*N.z # time varying force
  190. B1.apply_force(force, P1, B2, P2) #applying equal and opposite force on moving points
  191. assert B1.loads == [(P1, force)]
  192. assert B2.loads == [(P2, -force)]
  193. g1 = B1.mass*g*N.y
  194. g2 = B2.mass*g*N.y
  195. B1.apply_force(g1) #applying gravity on B1 masscenter
  196. B2.apply_force(g2) #applying gravity on B2 masscenter
  197. assert B1.loads == [(P1,force), (B1.masscenter, g1)]
  198. assert B2.loads == [(P2, -force), (B2.masscenter, g2)]
  199. force2 = x*N.x
  200. B1.apply_force(force2, reaction_body=B2) #Applying time varying force on masscenter
  201. assert B1.loads == [(P1, force), (B1.masscenter, force2+g1)]
  202. assert B2.loads == [(P2, -force), (B2.masscenter, -force2+g2)]
  203. def test_apply_torque():
  204. t = symbols('t')
  205. q = dynamicsymbols('q')
  206. with warns_deprecated_sympy():
  207. B1 = Body('B1')
  208. B2 = Body('B2')
  209. N = ReferenceFrame('N')
  210. torque = t*q*N.x
  211. B1.apply_torque(torque, B2) #Applying equal and opposite torque
  212. assert B1.loads == [(B1.frame, torque)]
  213. assert B2.loads == [(B2.frame, -torque)]
  214. torque2 = t*N.y
  215. B1.apply_torque(torque2)
  216. assert B1.loads == [(B1.frame, torque+torque2)]
  217. def test_clear_load():
  218. a = symbols('a')
  219. P = Point('P')
  220. with warns_deprecated_sympy():
  221. B = Body('B')
  222. force = a*B.z
  223. B.apply_force(force, P)
  224. assert B.loads == [(P, force)]
  225. B.clear_loads()
  226. assert B.loads == []
  227. def test_remove_load():
  228. P1 = Point('P1')
  229. P2 = Point('P2')
  230. with warns_deprecated_sympy():
  231. B = Body('B')
  232. f1 = B.x
  233. f2 = B.y
  234. B.apply_force(f1, P1)
  235. B.apply_force(f2, P2)
  236. assert B.loads == [(P1, f1), (P2, f2)]
  237. B.remove_load(P2)
  238. assert B.loads == [(P1, f1)]
  239. B.apply_torque(f1.cross(f2))
  240. assert B.loads == [(P1, f1), (B.frame, f1.cross(f2))]
  241. B.remove_load()
  242. assert B.loads == [(P1, f1)]
  243. def test_apply_loads_on_multi_degree_freedom_holonomic_system():
  244. """Example based on: https://pydy.readthedocs.io/en/latest/examples/multidof-holonomic.html"""
  245. with warns_deprecated_sympy():
  246. W = Body('W') #Wall
  247. B = Body('B') #Block
  248. P = Body('P') #Pendulum
  249. b = Body('b') #bob
  250. q1, q2 = dynamicsymbols('q1 q2') #generalized coordinates
  251. k, c, g, kT = symbols('k c g kT') #constants
  252. F, T = dynamicsymbols('F T') #Specified forces
  253. #Applying forces
  254. B.apply_force(F*W.x)
  255. W.apply_force(k*q1*W.x, reaction_body=B) #Spring force
  256. W.apply_force(c*q1.diff()*W.x, reaction_body=B) #dampner
  257. P.apply_force(P.mass*g*W.y)
  258. b.apply_force(b.mass*g*W.y)
  259. #Applying torques
  260. P.apply_torque(kT*q2*W.z, reaction_body=b)
  261. P.apply_torque(T*W.z)
  262. assert B.loads == [(B.masscenter, (F - k*q1 - c*q1.diff())*W.x)]
  263. assert P.loads == [(P.masscenter, P.mass*g*W.y), (P.frame, (T + kT*q2)*W.z)]
  264. assert b.loads == [(b.masscenter, b.mass*g*W.y), (b.frame, -kT*q2*W.z)]
  265. assert W.loads == [(W.masscenter, (c*q1.diff() + k*q1)*W.x)]
  266. def test_parallel_axis():
  267. N = ReferenceFrame('N')
  268. m, Ix, Iy, Iz, a, b = symbols('m, I_x, I_y, I_z, a, b')
  269. Io = inertia(N, Ix, Iy, Iz)
  270. # Test RigidBody
  271. o = Point('o')
  272. p = o.locatenew('p', a * N.x + b * N.y)
  273. with warns_deprecated_sympy():
  274. R = Body('R', masscenter=o, frame=N, mass=m, central_inertia=Io)
  275. Ip = R.parallel_axis(p)
  276. Ip_expected = inertia(N, Ix + m * b**2, Iy + m * a**2,
  277. Iz + m * (a**2 + b**2), ixy=-m * a * b)
  278. assert Ip == Ip_expected
  279. # Reference frame from which the parallel axis is viewed should not matter
  280. A = ReferenceFrame('A')
  281. A.orient_axis(N, N.z, 1)
  282. assert simplify(
  283. (R.parallel_axis(p, A) - Ip_expected).to_matrix(A)) == zeros(3, 3)
  284. # Test Particle
  285. o = Point('o')
  286. p = o.locatenew('p', a * N.x + b * N.y)
  287. with warns_deprecated_sympy():
  288. P = Body('P', masscenter=o, mass=m, frame=N)
  289. Ip = P.parallel_axis(p, N)
  290. Ip_expected = inertia(N, m * b ** 2, m * a ** 2, m * (a ** 2 + b ** 2),
  291. ixy=-m * a * b)
  292. assert not P.is_rigidbody
  293. assert Ip == Ip_expected