2import pytest
4from cymsim.config import Material, Numerics, Profile, Imperfections
5from cymsim.geometry import CymbalProfile
6from cymsim.mesh import disk_mesh
7from cymsim.shell import ShellModel
11def shell():
12 prof = CymbalProfile(Profile(), Imperfections(n_dents=30))
13 mesh = disk_mesh(Profile.radius, 0.025)
14 return ShellModel(mesh, prof, Material(), Numerics(dtype="float64"))
17def test_mass_matches_integrated_density(shell):
18 rho = shell.material.rho
19 expected = rho * (shell.h_e * shell.mesh.area).sum()
20 assert np.isclose(shell.mass, expected)
21 assert (shell.M > 0).all()
24def test_element_checks(shell):
25 c = shell.element_checks()
26 assert c.rigid_tilt_residual < 1e-12
27 assert c.constant_curvature_error < 1e-8
28 assert c.zero_energy_modes_min == 6 and c.zero_energy_modes_max == 6
29 assert c.rigid_mode_ratio < 1e-10
32def test_rigid_modes_are_M_orthonormal_and_force_free(shell):
33 import jax
34 import jax.numpy as jnp
35 G = np.einsum("knd,nd,lnd->kl", shell.Rm, shell.M, shell.Rm)
36 assert np.allclose(G, np.eye(6), atol=1e-10)
37 # the tangent stiffness at rest annihilates every rigid mode
38 q0 = jnp.zeros(shell.M.shape)
39 K = lambda d: np.asarray(jax.jvp(shell.internal_force, (q0,), (jnp.asarray(d),))[1])
40 rng = np.random.default_rng(0)
41 for k in range(6):
42 rnd = rng.standard_normal(shell.M.shape)
43 rnd *= np.linalg.norm(shell.Rm[k]) / np.linalg.norm(rnd)
44 assert np.abs(K(shell.Rm[k])).max() < 1e-12 * np.abs(K(rnd)).max()
45 # a finite tilt is force-free only to first order: the residual scales with amplitude squared
46 r1 = np.abs(np.asarray(shell.internal_force(shell.Rm_j[4] * 1e-3))).max()
47 r2 = np.abs(np.asarray(shell.internal_force(shell.Rm_j[4] * 2e-3))).max()
48 assert 3.5 < r2 / r1 < 4.5
51def test_projection_removes_rigid_motion(shell):
52 rng = np.random.default_rng(1)
53 v = rng.standard_normal(shell.M.shape)
54 p = np.asarray(shell.project_rigid(v))
55 assert np.allclose(np.einsum("knd,nd->k", shell.Rm, shell.M * p), 0, atol=1e-9)
56 assert np.allclose(np.asarray(shell.project_rigid(p)), p)
59def test_internal_force_is_minus_energy_gradient(shell):
60 rng = np.random.default_rng(2)
61 q = 1e-4 * rng.standard_normal(shell.M.shape)
62 d = rng.standard_normal(shell.M.shape)
63 eps = 1e-6
64 dE = (float(shell.total_energy(q + eps * d)) - float(shell.total_energy(q - eps * d))) / (2 * eps)
65 assert np.isclose(-(np.asarray(shell.internal_force(q)) * d).sum(), dE, rtol=1e-5)
68def test_damping_is_dissipative(shell):
69 rng = np.random.default_rng(3)
70 vw = rng.standard_normal(shell.mesh.n_nodes)
71 fd = np.asarray(shell.damping_force_w(vw))
72 assert (vw * fd).sum() > 0 # the force -fd does negative work on vw
73 assert np.isclose(fd.sum(), 0.0, atol=1e-9 * np.abs(fd).sum()) # a uniform velocity is undamped
76def test_stability_bound_is_finite_and_sane(shell):
77 assert 1e4 < shell.omega_max < 1e8
78 assert shell.dt_stable == 2.0 / shell.omega_max