diff --git a/elastica/rigidbody/mesh_rigid_body.py b/elastica/rigidbody/mesh_rigid_body.py index 509218f95..e3828fbcb 100644 --- a/elastica/rigidbody/mesh_rigid_body.py +++ b/elastica/rigidbody/mesh_rigid_body.py @@ -34,9 +34,17 @@ def __init__( self.density = density self.volume = volume self.mass = np.array([self.volume * self.density]) + self.mass_second_moment_of_inertia = mass_second_moment_of_inertia.reshape( + MaxDimension.value(), MaxDimension.value(), 1 + ) + + self.inv_mass_second_moment_of_inertia = np.linalg.inv( + mass_second_moment_of_inertia + ).reshape(MaxDimension.value(), MaxDimension.value(), 1) normal = np.array([1.0, 0.0, 0.0]).reshape(3, 1) tangents = np.array([0.0, 0.0, 1.0]).reshape(3, 1) binormal = _batch_cross(tangents, normal) + # initialize material frame to be the same as lab frame self.director_collection = np.zeros( (MaxDimension.value(), MaxDimension.value(), 1) @@ -48,34 +56,39 @@ def __init__( self.n_faces = mesh.faces.shape[-1] self.face_centers = np.array(mesh.face_centers.copy()) self.face_normals = np.array(mesh.face_normals.copy()) - self.distance_to_face_centers = _batch_norm( + + # since material frame is the same as lab frame initially, no need to convert this from lab to material + self.distance_to_face_centers_from_center_of_mass = _batch_norm( self.face_centers - center_of_mass.reshape(3, 1) ) - self.direction_to_face_centers = ( + self.direction_to_face_centers_from_center_of_mass_in_material_frame = ( self.face_centers - center_of_mass.reshape(3, 1) - ) / self.distance_to_face_centers - self.distance_to_faces = np.zeros((MaxDimension.value(), self.n_faces)) - self.direction_to_faces = np.zeros( + ) / self.distance_to_face_centers_from_center_of_mass + self.distance_to_faces_from_center_of_mass = np.zeros( + (MaxDimension.value(), self.n_faces) + ) + self.direction_to_faces_from_center_of_mass_in_material_frame = np.zeros( (MaxDimension.value(), MaxDimension.value(), self.n_faces) ) for i in range(MaxDimension.value()): for k in range(self.n_faces): - self.distance_to_faces[i, k] = np.linalg.norm( + self.distance_to_faces_from_center_of_mass[i, k] = np.linalg.norm( self.faces[:, i, k] - center_of_mass ) - self.direction_to_faces[:, i, k] = ( + self.direction_to_faces_from_center_of_mass_in_material_frame[ + :, i, k + ] = ( self.faces[:, i, k] - center_of_mass - ) / self.distance_to_faces[i, k] + ) / self.distance_to_faces_from_center_of_mass[ + i, k + ] - self.mass_second_moment_of_inertia = mass_second_moment_of_inertia.reshape( - MaxDimension.value(), MaxDimension.value(), 1 + self.face_normals_in_material_frame = np.zeros( + (MaxDimension.value(), self.n_faces) ) + self.face_normals_in_material_frame = self.face_normals.copy() - self.inv_mass_second_moment_of_inertia = np.linalg.inv( - mass_second_moment_of_inertia - ).reshape(MaxDimension.value(), MaxDimension.value(), 1) - - # position is at the center + # position is at the center of mass self.position_collection = np.zeros((MaxDimension.value(), 1)) self.position_collection[:, 0] = center_of_mass @@ -83,16 +96,6 @@ def __init__( self.omega_collection = np.zeros((MaxDimension.value(), 1)) self.acceleration_collection = np.zeros((MaxDimension.value(), 1)) self.alpha_collection = np.zeros((MaxDimension.value(), 1)) - self.face_normals_lagrangian = np.zeros((MaxDimension.value(), self.n_faces)) - - # self.face_normals_lagrangian = self.director_collection[...,0].T@self.face_normals - # for i in range(3): - # for j in range(3): - # for k in range(self.n_faces): - # self.face_normals_lagrangian[i, k] += ( - # self.director_collection[j, i, 0] * self.face_normals[j, k] - # ) - self.face_normals_lagrangian = self.face_normals.copy() self.external_forces = np.zeros((MaxDimension.value())).reshape( MaxDimension.value(), 1 @@ -106,22 +109,22 @@ def update_faces(self): self.director_collection, self.faces, self.position_collection, - self.distance_to_faces, - self.direction_to_faces, + self.distance_to_faces_from_center_of_mass, + self.direction_to_faces_from_center_of_mass_in_material_frame, self.n_faces, ) _update_face_centers( self.director_collection, self.face_centers, self.position_collection, - self.distance_to_face_centers, - self.direction_to_face_centers, + self.distance_to_face_centers_from_center_of_mass, + self.direction_to_face_centers_from_center_of_mass_in_material_frame, self.n_faces, ) _update_face_normals( self.director_collection, self.face_normals, - self.face_normals_lagrangian, + self.face_normals_in_material_frame, self.n_faces, ) @@ -131,8 +134,8 @@ def _update_face_centers( director_collection, face_centers, center_of_mass, - distance_to_face_centers, - direction_to_face_centers, + distance_to_face_centers_from_center_of_mass, + direction_to_face_centers_from_center_of_mass_in_material_frame, n_faces, ): face_centers[:] = np.zeros((3, n_faces)) @@ -141,9 +144,11 @@ def _update_face_centers( face_centers[i, k] += center_of_mass[i, 0] for j in range(3): face_centers[i, k] += ( - distance_to_face_centers[i] + distance_to_face_centers_from_center_of_mass[i] * director_collection[i, j, 0] - * direction_to_face_centers[j, k] + * direction_to_face_centers_from_center_of_mass_in_material_frame[ + j, k + ] ) @@ -152,8 +157,8 @@ def _update_faces( director_collection, faces, center_of_mass, - distance_to_faces, - direction_to_faces, + distance_to_faces_from_center_of_mass, + direction_to_faces_from_center_of_mass_in_material_frame, n_faces, ): faces[:] = np.zeros((3, 3, n_faces)) # dim,vertices,faces @@ -163,9 +168,11 @@ def _update_faces( faces[i, m, k] += center_of_mass[i, 0] for j in range(3): faces[i, m, k] += ( - distance_to_faces[m, k] + distance_to_faces_from_center_of_mass[m, k] * director_collection[i, j, 0] - * direction_to_faces[j, m, k] + * direction_to_faces_from_center_of_mass_in_material_frame[ + j, m, k + ] ) @@ -173,7 +180,7 @@ def _update_faces( def _update_face_normals( director_collection, face_normals, - face_normals_lagrangian, + face_normals_in_material_frame, n_faces, ): @@ -182,5 +189,5 @@ def _update_face_normals( for j in range(3): for k in range(n_faces): face_normals[i, k] += ( - director_collection[i, j, 0] * face_normals_lagrangian[j, k] + director_collection[i, j, 0] * face_normals_in_material_frame[j, k] ) diff --git a/tests/test_rigid_body/test_mesh_rigid_body.py b/tests/test_rigid_body/test_mesh_rigid_body.py index cf6e3948a..f7a54cb3c 100644 --- a/tests/test_rigid_body/test_mesh_rigid_body.py +++ b/tests/test_rigid_body/test_mesh_rigid_body.py @@ -4,165 +4,175 @@ from elastica.utils import Tolerance from elastica.rigidbody import MeshRigidBody +from elastica import Mesh from elastica._linalg import ( - _batch_cross, + _batch_norm, ) -class MockMesh: - def __init__(self, n_faces, faces, face_centers, face_normals): - self.n_faces = n_faces - self.faces = faces - self.face_centers = face_centers - self.face_normals = face_normals +def initialize_cube_rigid_body(): + """ + This function is to initialize the cube rigid body from the cube.stl. + """ + cube_mesh = Mesh(r"tests/cube.stl") + center_of_mass = np.array([0.0, 0.0, 0.0]) + base_length = 2 + volume = base_length ** 3 + density = 1.0 + mass = density * volume + # Mass second moment of inertia for cube + mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) + np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) + cube_mesh_rigid_body = MeshRigidBody( + cube_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume + ) + return ( + cube_mesh_rigid_body, + cube_mesh, + center_of_mass, + mass, + mass_second_moment_of_inertia, + ) -# tests Initialisation of mesh rigid body +# tests Initialization of mesh rigid body def test_MeshRigidBody_initialization(): """ This test case is for testing initialization of mesh rigid body and it checks the validity of the members of MeshRigidBody class. - - Returns - ------- - """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] + ( + cube_mesh_rigid_body, + cube_mesh, + correct_center_of_mass, + correct_mass, + correct_mass_second_moment_of_inertia, + ) = initialize_cube_rigid_body() + + # checking mesh rigid body center + assert_allclose( + cube_mesh_rigid_body.position_collection[..., -1], + correct_center_of_mass, + atol=Tolerance.atol(), ) - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] + # checking velocities, omegas, accelerations and alphas + assert_allclose( + cube_mesh_rigid_body.velocity_collection, + np.zeros((3, 1)), + atol=Tolerance.atol(), ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] + assert_allclose( + cube_mesh_rigid_body.omega_collection, + np.zeros((3, 1)), + atol=Tolerance.atol(), + ) + assert_allclose( + cube_mesh_rigid_body.acceleration_collection, + np.zeros((3, 1)), + atol=Tolerance.atol(), + ) + assert_allclose( + cube_mesh_rigid_body.alpha_collection, + np.zeros((3, 1)), + atol=Tolerance.atol(), ) - n_faces = 12 - # setting up test params - test_mesh = MockMesh(n_faces, faces, face_centers, face_normals) - center_of_mass = np.array([0.0, 0.0, 0.0]) - direction = np.array([0.0, 0.0, 1.0]).reshape(3, 1) - normal = np.array([1.0, 0.0, 0.0]).reshape(3, 1) - binormal = np.array([0.0, 1.0, 0.0]).reshape(3, 1) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - # Inverse mass second of inertia - inv_mass_second_moment_of_inertia = np.linalg.inv(mass_second_moment_of_inertia) + # check mass and density initalization + cube_base_length = 2 + cube_volume = cube_base_length ** 3 + correct_density = 1.0 + assert_allclose( + cube_mesh_rigid_body.density, correct_density, atol=Tolerance.atol() + ) + assert_allclose(cube_mesh_rigid_body.volume, cube_volume, atol=Tolerance.atol()) + assert_allclose(cube_mesh_rigid_body.mass, correct_mass, atol=Tolerance.atol()) - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume + # check mass second moment of inertia initalization + correct_inv_mass_second_moment_of_inertia = np.linalg.inv( + correct_mass_second_moment_of_inertia ) - # checking origin and length of rod assert_allclose( - test_mesh_rigid_body.position_collection[..., -1], - center_of_mass, + cube_mesh_rigid_body.inv_mass_second_moment_of_inertia[..., -1], + correct_inv_mass_second_moment_of_inertia, atol=Tolerance.atol(), ) - - # element lengths are equal for all rod. - # checking velocities, omegas and rest strains - # density and mass assert_allclose( - test_mesh_rigid_body.velocity_collection, - np.zeros((3, 1)), + cube_mesh_rigid_body.mass_second_moment_of_inertia[..., -1], + correct_mass_second_moment_of_inertia, atol=Tolerance.atol(), ) + # check director initalization correct_director_collection = np.zeros((3, 3, 1)) - correct_director_collection[0] = normal - correct_director_collection[1] = binormal - correct_director_collection[2] = direction + correct_director_collection[0] = np.array([1.0, 0.0, 0.0]).reshape(3, 1) + correct_director_collection[1] = np.array([0.0, 1.0, 0.0]).reshape(3, 1) + correct_director_collection[2] = np.array([0.0, 0.0, 1.0]).reshape(3, 1) assert_allclose( - test_mesh_rigid_body.director_collection, + cube_mesh_rigid_body.director_collection, correct_director_collection, atol=Tolerance.atol(), ) + # check faces, n_faces, face_centers, face_normals + correct_faces = cube_mesh.faces + correct_n_faces = cube_mesh.faces.shape[-1] + correct_face_centers = cube_mesh.face_centers + correct_face_normals = cube_mesh.face_normals + assert_allclose(cube_mesh_rigid_body.faces, correct_faces, atol=Tolerance.atol()) + assert_allclose( + cube_mesh_rigid_body.n_faces, correct_n_faces, atol=Tolerance.atol() + ) + assert_allclose( + cube_mesh_rigid_body.face_centers, correct_face_centers, atol=Tolerance.atol() + ) assert_allclose( - test_mesh_rigid_body.omega_collection, np.zeros((3, 1)), atol=Tolerance.atol() + cube_mesh_rigid_body.face_normals, correct_face_normals, atol=Tolerance.atol() ) - assert_allclose(test_mesh_rigid_body.density, density, atol=Tolerance.atol()) + # check distance to faces, direction to faces, distance to face center, direction to faces centers, face normals in material frame + correct_distance_to_face_centers_from_center_of_mass = _batch_norm( + correct_face_centers - correct_center_of_mass.reshape(3, 1) + ) + correct_direction_to_face_centers_from_center_of_mass_in_material_frame = ( + correct_face_centers - correct_center_of_mass.reshape(3, 1) + ) / correct_distance_to_face_centers_from_center_of_mass - # Check mass at each node. Note that, node masses is - # half of element mass at the first and last node. - assert_allclose(test_mesh_rigid_body.mass, mass, atol=Tolerance.atol()) + correct_distance_to_faces_from_center_of_mass = np.zeros((3, correct_n_faces)) + correct_direction_to_faces_from_center_of_mass_in_material_frame = np.zeros( + (3, 3, correct_n_faces) + ) + for i in range(3): + for k in range(correct_n_faces): + correct_distance_to_faces_from_center_of_mass[i, k] = np.linalg.norm( + correct_faces[:, i, k] - correct_center_of_mass + ) + correct_direction_to_faces_from_center_of_mass_in_material_frame[ + :, i, k + ] = ( + correct_faces[:, i, k] - correct_center_of_mass + ) / correct_distance_to_faces_from_center_of_mass[ + i, k + ] - # checking directors, rest length - # and shear, bend matrices and moment of inertia assert_allclose( - test_mesh_rigid_body.inv_mass_second_moment_of_inertia[..., -1], - inv_mass_second_moment_of_inertia, + cube_mesh_rigid_body.distance_to_face_centers_from_center_of_mass, + correct_distance_to_face_centers_from_center_of_mass, + atol=Tolerance.atol(), + ) + assert_allclose( + cube_mesh_rigid_body.direction_to_face_centers_from_center_of_mass_in_material_frame, + correct_direction_to_face_centers_from_center_of_mass_in_material_frame, + atol=Tolerance.atol(), + ) + assert_allclose( + cube_mesh_rigid_body.distance_to_faces_from_center_of_mass, + correct_distance_to_faces_from_center_of_mass, + atol=Tolerance.atol(), + ) + assert_allclose( + cube_mesh_rigid_body.direction_to_faces_from_center_of_mass_in_material_frame, + correct_direction_to_faces_from_center_of_mass_in_material_frame, atol=Tolerance.atol(), ) @@ -170,1122 +180,168 @@ def test_MeshRigidBody_initialization(): def test_mesh_rigid_body_update_accelerations(): """ This test is testing the update acceleration method of MeshRigidBody class. - - Returns - ------- - """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - n_faces = 12 - - test_mesh = MockMesh(n_faces, faces, face_centers, face_normals) - center_of_mass = np.zeros( - 3, - ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume + ( + cube_mesh_rigid_body, + cube_mesh, + correct_center_of_mass, + correct_mass, + correct_mass_second_moment_of_inertia, + ) = initialize_cube_rigid_body() # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume - ) - - inv_mass_second_moment_of_inertia = ( - test_mesh_rigid_body.inv_mass_second_moment_of_inertia.reshape(3, 3) + inv_mass_second_moment_of_inertia = np.linalg.inv( + correct_mass_second_moment_of_inertia ) - external_forces = np.random.randn(3).reshape(3, 1) external_torques = np.random.randn(3).reshape(3, 1) - correct_acceleration = external_forces / mass + correct_acceleration = external_forces / correct_mass correct_alpha = inv_mass_second_moment_of_inertia @ external_torques.reshape(3) correct_alpha = correct_alpha.reshape(3, 1) - test_mesh_rigid_body.external_forces[:] = external_forces - test_mesh_rigid_body.external_torques[:] = external_torques + cube_mesh_rigid_body.external_forces[:] = external_forces + cube_mesh_rigid_body.external_torques[:] = external_torques - test_mesh_rigid_body.update_accelerations(time=0) + cube_mesh_rigid_body.update_accelerations(time=0) assert_allclose( correct_acceleration, - test_mesh_rigid_body.acceleration_collection, + cube_mesh_rigid_body.acceleration_collection, atol=Tolerance.atol(), ) assert_allclose( - correct_alpha, test_mesh_rigid_body.alpha_collection, atol=Tolerance.atol() + correct_alpha, cube_mesh_rigid_body.alpha_collection, atol=Tolerance.atol() ) def test_compute_position_center_of_mass(): """ This test is testing compute position center of mass method of MeshRigidBody class. - - Returns - ------- - """ - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - n_faces = 12 - - test_mesh = MockMesh(n_faces, faces, face_centers, face_normals) - center_of_mass = np.zeros( - 3, - ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume + ( + cube_mesh_rigid_body, + cube_mesh, + correct_center_of_mass, + correct_mass, + correct_mass_second_moment_of_inertia, + ) = initialize_cube_rigid_body() + assert_allclose( + correct_center_of_mass, + cube_mesh_rigid_body.compute_position_center_of_mass(), + atol=Tolerance.atol(), ) - correct_position = center_of_mass - - test_position = test_mesh_rigid_body.compute_position_center_of_mass() - - assert_allclose(correct_position, test_position, atol=Tolerance.atol()) - def test_compute_translational_energy(): """ This test is testing compute translational energy function. - Returns - ------- - """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - n_faces = 12 - - test_mesh = MockMesh(n_faces, faces, face_centers, face_normals) - center_of_mass = np.zeros( - 3, - ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume - ) - + ( + cube_mesh_rigid_body, + cube_mesh, + correct_center_of_mass, + correct_mass, + correct_mass_second_moment_of_inertia, + ) = initialize_cube_rigid_body() speed = np.random.randn() - test_mesh_rigid_body.velocity_collection[2] = speed - - correct_energy = 0.5 * mass * speed ** 2 - test_energy = test_mesh_rigid_body.compute_translational_energy() - - assert_allclose(correct_energy, test_energy, atol=Tolerance.atol()) - - -def test_update_faces_translation(): - """ - This test is testing compute position center of mass method of MeshRigidBody class. - - Returns - ------- - - """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - n_faces = 12 - - test_mesh = MockMesh(n_faces, faces, face_centers, face_normals) - center_of_mass = np.zeros( - 3, - ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume - ) - - test_mesh_rigid_body.position_collection[:, 0] = np.ones( - 3, - ) - - correct_face_centers = face_centers + 1.0 * np.ones_like(face_centers) - correct_faces = faces + 1.0 * np.ones_like(faces) - correct_face_normals = face_normals - test_mesh_rigid_body.update_faces() - test_face_centers = test_mesh_rigid_body.face_centers - test_face_normals = test_mesh_rigid_body.face_normals - test_faces = test_mesh_rigid_body.faces - - assert_allclose(correct_face_centers, test_face_centers, atol=Tolerance.atol()) - assert_allclose(correct_faces, test_faces, atol=Tolerance.atol()) - assert_allclose(correct_face_normals, test_face_normals, atol=Tolerance.atol()) - - -def test_update_faces_rotation(): - """ - This test is testing compute position center of mass method of MeshRigidBody class. - - Returns - ------- - - """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - n_faces = 12 - - test_mesh = MockMesh(n_faces, faces, face_centers, face_normals) - center_of_mass = np.zeros( - 3, - ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume - ) - - normal = np.array([0.0, 0.0, 1.0]).reshape(3, 1) - tangent = np.array([-1.0, 0.0, 0.0]).reshape(3, 1) - binormal = _batch_cross(tangent, normal) - - rotation_matrix = np.zeros((3, 3)) - rotation_matrix[0, ...] = normal.reshape( - 3, - ) - rotation_matrix[1, ...] = binormal.reshape( - 3, - ) - rotation_matrix[2, ...] = tangent.reshape( - 3, - ) - correct_faces = np.zeros_like(faces) - - correct_face_centers = rotation_matrix @ face_centers - for i in range(3): - correct_faces[:, i, :] = rotation_matrix @ faces[:, i, :] - correct_face_normals = rotation_matrix @ face_normals - - test_mesh_rigid_body.director_collection[0, ...] = normal - test_mesh_rigid_body.director_collection[1, ...] = binormal - test_mesh_rigid_body.director_collection[2, ...] = tangent - test_mesh_rigid_body.update_faces() - - test_face_centers = test_mesh_rigid_body.face_centers - test_face_normals = test_mesh_rigid_body.face_normals - test_faces = test_mesh_rigid_body.faces - - assert_allclose(correct_face_normals, test_face_normals, atol=Tolerance.atol()) - assert_allclose(correct_face_centers, test_face_centers, atol=Tolerance.atol()) - assert_allclose(correct_faces, test_faces, atol=Tolerance.atol()) - - -def test_update_faces_rotation_and_translation(): - """ - This test is testing compute position center of mass method of MeshRigidBody class. - - Returns - ------- - - """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - n_faces = 12 - - test_mesh = MockMesh(n_faces, faces, face_centers, face_normals) - center_of_mass = np.zeros( - 3, - ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume - ) - - normal = np.array([1.0 / np.sqrt(2), -1.0 / np.sqrt(2), 0.0]).reshape(3, 1) - tangent = np.array([1.0 / np.sqrt(2), 1.0 / np.sqrt(2), 0.0]).reshape(3, 1) - binormal = _batch_cross(tangent, normal) - - rotation_matrix = np.zeros((3, 3)) - rotation_matrix[0, ...] = normal.reshape( - 3, - ) - rotation_matrix[1, ...] = binormal.reshape( - 3, - ) - rotation_matrix[2, ...] = tangent.reshape( - 3, - ) - - translation = np.ones_like(face_centers) - correct_faces = np.zeros_like(faces) - for i in range(3): - correct_faces[:, i, :] = translation + rotation_matrix @ faces[:, i, :] - correct_face_centers = translation + rotation_matrix @ face_centers - correct_face_normals = rotation_matrix @ face_normals - - test_mesh_rigid_body.director_collection[0, ...] = normal - test_mesh_rigid_body.director_collection[1, ...] = binormal - test_mesh_rigid_body.director_collection[2, ...] = tangent - test_mesh_rigid_body.position_collection[:, 0] = np.ones( - 3, + cube_mesh_rigid_body.velocity_collection[2] = speed + assert_allclose( + 0.5 * correct_mass * speed ** 2, + cube_mesh_rigid_body.compute_translational_energy(), + atol=Tolerance.atol(), ) - print(test_mesh_rigid_body.face_normals) - test_mesh_rigid_body.update_faces() - test_face_centers = test_mesh_rigid_body.face_centers.copy() - test_face_normals = test_mesh_rigid_body.face_normals.copy() - print(test_face_normals) - test_faces = test_mesh_rigid_body.faces.copy() - assert_allclose(correct_face_normals, test_face_normals, atol=Tolerance.atol()) - assert_allclose(correct_face_centers, test_face_centers, atol=Tolerance.atol()) - assert_allclose(correct_faces, test_faces, atol=Tolerance.atol()) - - -def test_update_faces_translation_with_mesh_initializer(): +def test_update_faces_translation(): """ - This test is testing compute position center of mass method of MeshRigidBody class. - - Returns - ------- - + This test is testing update_faces method of MeshRigidBody class (translation only). """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - from elastica import Mesh - - test_mesh = Mesh(r"tests/cube.stl") - center_of_mass = np.zeros( + ( + cube_mesh_rigid_body, + cube_mesh, + correct_center_of_mass, + correct_mass, + correct_mass_second_moment_of_inertia, + ) = initialize_cube_rigid_body() + + cube_mesh_rigid_body.position_collection[:, 0] = np.ones( 3, ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume - ) - translation = np.ones_like(face_centers) - correct_faces = np.zeros_like(faces) - for i in range(3): - correct_faces[:, i, :] = translation + faces[:, i, :] - correct_face_centers = translation + face_centers - correct_face_normals = face_normals - - test_mesh_rigid_body.position_collection[:, 0] = np.ones( - 3, + correct_face_centers = cube_mesh.face_centers + 1.0 * np.ones_like( + cube_mesh.face_centers ) - test_mesh_rigid_body.update_faces() + correct_faces = cube_mesh.faces + 1.0 * np.ones_like(cube_mesh.faces) + correct_face_normals = cube_mesh.face_normals + cube_mesh_rigid_body.update_faces() - test_face_centers = test_mesh_rigid_body.face_centers - test_face_normals = test_mesh_rigid_body.face_normals - test_faces = test_mesh_rigid_body.faces - - assert_allclose( - correct_face_normals, - test_face_normals, - atol=Tolerance.atol(), - rtol=Tolerance.rtol(), - ) assert_allclose( - correct_face_centers, - test_face_centers, - atol=Tolerance.atol(), - rtol=Tolerance.rtol(), + correct_face_normals, cube_mesh_rigid_body.face_normals, atol=Tolerance.atol() ) assert_allclose( - correct_faces, test_faces, atol=Tolerance.atol(), rtol=Tolerance.rtol() + correct_face_centers, cube_mesh_rigid_body.face_centers, atol=Tolerance.atol() ) + assert_allclose(correct_faces, cube_mesh_rigid_body.faces, atol=Tolerance.atol()) -def test_update_faces_rotation_with_mesh_initializer(): +def test_update_faces_rotation(): """ - This test is testing compute position center of mass method of MeshRigidBody class. - - Returns - ------- - + This test is testing update_faces method of MeshRigidBody class (rotation only). """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - from elastica import Mesh - - test_mesh = Mesh(r"tests/cube.stl") - center_of_mass = np.zeros( - 3, - ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume - ) - - normal = np.array([0.0, 0.0, 1.0]).reshape(3, 1) - tangent = np.array([-1.0, 0.0, 0.0]).reshape(3, 1) - binormal = _batch_cross(tangent, normal) + ( + cube_mesh_rigid_body, + cube_mesh, + correct_center_of_mass, + correct_mass, + correct_mass_second_moment_of_inertia, + ) = initialize_cube_rigid_body() + normal = np.array([0.0, 0.0, 1.0]) + tangent = np.array([-1.0, 0.0, 0.0]) + binormal = np.cross(tangent, normal) rotation_matrix = np.zeros((3, 3)) - rotation_matrix[0, ...] = normal.reshape( - 3, - ) - rotation_matrix[1, ...] = binormal.reshape( - 3, - ) - rotation_matrix[2, ...] = tangent.reshape( - 3, - ) + rotation_matrix[0, ...] = normal + rotation_matrix[1, ...] = binormal + rotation_matrix[2, ...] = tangent - correct_faces = np.zeros_like(faces) + correct_faces = np.zeros_like(cube_mesh.faces) for i in range(3): - correct_faces[:, i, :] = rotation_matrix @ faces[:, i, :] - correct_face_centers = rotation_matrix @ test_mesh.face_centers - correct_face_normals = rotation_matrix @ test_mesh.face_normals + correct_faces[:, i, :] = rotation_matrix @ cube_mesh.faces[:, i, :] + correct_face_centers = rotation_matrix @ cube_mesh.face_centers + correct_face_normals = rotation_matrix @ cube_mesh.face_normals - test_mesh_rigid_body.director_collection[0, ...] = normal - test_mesh_rigid_body.director_collection[1, ...] = binormal - test_mesh_rigid_body.director_collection[2, ...] = tangent - test_mesh_rigid_body.update_faces() + cube_mesh_rigid_body.director_collection[0, ...] = normal.reshape(3, 1) + cube_mesh_rigid_body.director_collection[1, ...] = binormal.reshape(3, 1) + cube_mesh_rigid_body.director_collection[2, ...] = tangent.reshape(3, 1) + cube_mesh_rigid_body.update_faces() - test_face_centers = test_mesh_rigid_body.face_centers - test_face_normals = test_mesh_rigid_body.face_normals - test_faces = test_mesh_rigid_body.faces - - assert_allclose( - correct_face_normals, - test_face_normals, - atol=Tolerance.atol(), - rtol=Tolerance.rtol(), - ) assert_allclose( - correct_face_centers, - test_face_centers, - atol=Tolerance.atol(), - rtol=Tolerance.rtol(), + correct_face_normals, cube_mesh_rigid_body.face_normals, atol=Tolerance.atol() ) assert_allclose( - correct_faces, test_faces, atol=Tolerance.atol(), rtol=Tolerance.rtol() + correct_face_centers, cube_mesh_rigid_body.face_centers, atol=Tolerance.atol() ) + assert_allclose(correct_faces, cube_mesh_rigid_body.faces, atol=Tolerance.atol()) -def test_update_faces_rotation_and_translation_with_mesh_initializer(): +def test_update_faces_rotation_and_translation(): """ - This test is testing compute position center of mass method of MeshRigidBody class. - - Returns - ------- - + This test is testing update_faces method of MeshRigidBody class (rotation and translation). """ - - faces = 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], - [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], - [-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, 1.0], - ], - ] - ) - - face_normals = np.array( - [ - [0, 0, -1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0], - [0, 0, 0, 0, 0, 0, 0, 0, 1.0, 1.0, -1.0, -1.0], - [-1.0, -1.0, 0, 0, 1.0, 1.0, 0, 0, 0, 0, 0, 0], - ] - ) - face_centers = np.array( - [ - [ - -0.333333, - 0.333333, - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - ], - [ - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - 1.0, - 1.0, - -1.0, - -1.0, - ], - [ - -1.0, - -1.0, - 0.333333, - -0.333333, - 1.0, - 1.0, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - -0.333333, - 0.333333, - ], - ] - ) - from elastica import Mesh - - test_mesh = Mesh(r"tests/cube.stl") - center_of_mass = np.zeros( - 3, - ) - base_length = 2 - volume = base_length ** 3 - density = np.random.uniform(1.0, 10) - mass = density * volume - # Mass second moment of inertia for cube - mass_second_moment_of_inertia = np.zeros((3, 3), np.float64) - np.fill_diagonal(mass_second_moment_of_inertia, (mass * base_length ** 2) / 6) - - test_mesh_rigid_body = MeshRigidBody( - test_mesh, center_of_mass, mass_second_moment_of_inertia, density, volume - ) - - normal = np.array([0.0, 0.0, 1.0]).reshape(3, 1) - tangent = np.array([-1.0, 0.0, 0.0]).reshape(3, 1) - binormal = _batch_cross(tangent, normal) + ( + cube_mesh_rigid_body, + cube_mesh, + correct_center_of_mass, + correct_mass, + correct_mass_second_moment_of_inertia, + ) = initialize_cube_rigid_body() + + normal = np.array([1.0 / np.sqrt(2), -1.0 / np.sqrt(2), 0.0]) + tangent = np.array([1.0 / np.sqrt(2), 1.0 / np.sqrt(2), 0.0]) + binormal = np.cross(tangent, normal) rotation_matrix = np.zeros((3, 3)) rotation_matrix[0, ...] = normal.reshape( @@ -1298,40 +354,30 @@ def test_update_faces_rotation_and_translation_with_mesh_initializer(): 3, ) - translation = np.ones_like(face_centers) - correct_faces = np.zeros_like(faces) + translation = np.ones_like(cube_mesh.face_centers) + correct_faces = np.zeros_like(cube_mesh.faces) for i in range(3): - correct_faces[:, i, :] = translation + rotation_matrix @ faces[:, i, :] - correct_face_centers = translation + rotation_matrix @ face_centers - correct_face_normals = rotation_matrix @ face_normals - - test_mesh_rigid_body.director_collection[0, ...] = normal - test_mesh_rigid_body.director_collection[1, ...] = binormal - test_mesh_rigid_body.director_collection[2, ...] = tangent - test_mesh_rigid_body.position_collection[:, 0] = np.ones( + correct_faces[:, i, :] = ( + translation + rotation_matrix @ cube_mesh.faces[:, i, :] + ) + correct_face_centers = translation + rotation_matrix @ cube_mesh.face_centers + correct_face_normals = rotation_matrix @ cube_mesh.face_normals + + cube_mesh_rigid_body.director_collection[0, ...] = normal.reshape(3, 1) + cube_mesh_rigid_body.director_collection[1, ...] = binormal.reshape(3, 1) + cube_mesh_rigid_body.director_collection[2, ...] = tangent.reshape(3, 1) + cube_mesh_rigid_body.position_collection[:, 0] = np.ones( 3, ) - test_mesh_rigid_body.update_faces() - - test_face_centers = test_mesh_rigid_body.face_centers - test_face_normals = test_mesh_rigid_body.face_normals - test_faces = test_mesh_rigid_body.faces + cube_mesh_rigid_body.update_faces() assert_allclose( - correct_face_normals, - test_face_normals, - atol=Tolerance.atol(), - rtol=Tolerance.rtol(), - ) - assert_allclose( - correct_face_centers, - test_face_centers, - atol=Tolerance.atol(), - rtol=Tolerance.rtol(), + correct_face_normals, cube_mesh_rigid_body.face_normals, atol=Tolerance.atol() ) assert_allclose( - correct_faces, test_faces, atol=Tolerance.atol(), rtol=Tolerance.rtol() + correct_face_centers, cube_mesh_rigid_body.face_centers, atol=Tolerance.atol() ) + assert_allclose(correct_faces, cube_mesh_rigid_body.faces, atol=Tolerance.atol()) if __name__ == "__main__":