1import jax
2import jax.numpy as jnp
3import numpy as np
5from cymsim.elements import dkt_stiffness, element_energy_parts, membrane_strains
8def _triangle():
9 return jnp.array([0.0, 0.03, 0.01]), jnp.array([0.0, 0.005, 0.025])
12def test_dkt_rigid_modes_and_constant_curvature():
13 xe, ye = _triangle()
14 D, nu = 10.0, 0.3
15 Ke = np.asarray(dkt_stiffness(xe, ye, D, nu))
16 assert np.allclose(Ke, Ke.T)
17 x, y = np.asarray(xe), np.asarray(ye)
18 # rigid modes: translation, tilt about x (w = y, thx = 1), tilt about y (w = -x, thy = 1)
19 for U in (np.tile([1, 0, 0], 3), np.stack([y, np.ones(3), np.zeros(3)], 1).ravel(),
20 np.stack([-x, np.zeros(3), np.ones(3)], 1).ravel()):
21 assert np.abs(Ke @ U).max() < 1e-12 * np.abs(Ke).max()
22 # constant curvature w = x^2 / 2: energy D A / 2
23 A = 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]))
24 U = np.stack([0.5 * x**2, np.zeros(3), -x], 1).ravel()
25 assert np.isclose(0.5 * U @ Ke @ U, 0.5 * D * A, rtol=1e-9)
26 # twist w = x y: energy D (1 - nu) A
27 U = np.stack([x * y, x, -y], 1).ravel()
28 assert np.isclose(0.5 * U @ Ke @ U, D * (1 - nu) * A, rtol=1e-9)
29 eig = np.linalg.eigvalsh(Ke)
30 assert (eig[:3] < 1e-10 * eig[-1]).all() and eig[3] > 1e-6 * eig[-1]
33def test_membrane_strains_linear_and_nonlinear():
34 Gx, Gy = jnp.array([1.0, -1.0, 0.0]), jnp.array([0.0, 1.0, -1.0])
35 # nodal values arranged so that u_x = 2, v_y = 3, w_x = 0.1, w_y = 0.2 for these operators
36 qe = jnp.zeros((3, 5)).at[:, 0].set(jnp.array([2.0, 0.0, 0.0])).at[:, 1].set(jnp.array([0.0, 3.0, 0.0])) \
37 .at[:, 2].set(jnp.array([0.1, 0.0, -0.2]))
38 exx, eyy, gxy = membrane_strains(qe, Gx, Gy, 0.5, -0.5, 1.0)
39 assert np.isclose(exx, 2 + 0.5 * 0.01 + 0.5 * 0.1)
40 assert np.isclose(eyy, 3 + 0.5 * 0.04 - 0.5 * 0.2)
41 assert np.isclose(gxy, -3 + 0.1 * 0.2 + 0.5 * 0.2 - 0.5 * 0.1) # u_y + v_x = 0 - 3
42 exx0, eyy0, _ = membrane_strains(qe, Gx, Gy, 0.5, -0.5, 0.0)
43 assert np.isclose(exx0, 2 + 0.5 * 0.1) and np.isclose(eyy0, 3 - 0.5 * 0.2)
46def test_energy_hessian_is_symmetric_positive_semidefinite():
47 xe, ye = _triangle()
48 Ke = dkt_stiffness(xe, ye, 10.0, 0.3)
49 x, y = np.asarray(xe), np.asarray(ye)
50 A = 0.5 * ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0]))
51 Gx = jnp.array([y[1] - y[2], y[2] - y[0], y[0] - y[1]]) / (2 * A)
52 Gy = jnp.array([x[2] - x[1], x[0] - x[2], x[1] - x[0]]) / (2 * A)
53 Cm = jnp.array([1e6, 3e5, 3.5e5])
54 f = lambda q: sum(element_energy_parts(q, Gx, Gy, A, 0.1, -0.05, Cm, Ke, 1.0))
55 H = np.asarray(jax.hessian(f)(jnp.zeros((3, 5)))).reshape(15, 15)
56 assert np.allclose(H, H.T)
57 eig = np.linalg.eigvalsh(H)
58 assert eig[0] > -1e-9 * eig[-1]
59 assert (eig < 1e-9 * eig[-1]).sum() == 6 # 3 translations + 3 rotations