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
42c1bf8Time step from the global stability bound (power iteration), with numerics.stability_bound to keep the element boundJeremy Magland 78 assert shell.omega_global < shell.omega_max # the element bound is an upper bound
79 assert shell.dt_stable == 2.0 / shell.omega_global # the default numerics.stability_bound is "global"
82def test_global_bound_matches_the_dense_eigenproblem(shell):
83 """The power iteration finds the largest eigenvalue of M^-1 K, K assembled from the tangents at rest."""
84 N, L = shell.mesh.n_nodes, shell.M.size
85 K = np.zeros((L, L))
86 dof = (5 * shell.mesh.tri[:, :, None] + np.arange(5)).reshape(-1, 15)
87 for e in range(len(dof)):
88 K[np.ix_(dof[e], dof[e])] += shell.K_rest[e]
89 isq = 1 / np.sqrt(shell.M.ravel())
90 lam = np.linalg.eigvalsh(isq[:, None] * K * isq[None, :])[-1]
91 assert abs(shell.omega_global / np.sqrt(lam) - 1) < 1e-6
94def test_element_bound_on_request(shell):
95 sh = ShellModel(shell.mesh, CymbalProfile(Profile(), Imperfections(n_dents=30)), Material(),
96 Numerics(dtype="float64", stability_bound="element"))
97 assert sh.dt_stable == 2.0 / sh.omega_max