1"""Linear check of the assembly: the lowest natural frequencies of the linearized shell with the lumped and the
2consistent mass, from the same element operators the time stepping uses. For the flat annulus the free disk's
3values are 21.9 (2,0), 38.0 (0,0), 50.9 (3,0), 85.9 (1,0), 89.6 (4,0) Hz (Leissa), the 6 mm hole changing them
4by well under a percent. Usage: python check_modes.py <challenge.json> [spacing_m] [n_modes]"""
5import json
6import sys
8import numpy as np
9import scipy.sparse as sp
10import scipy.sparse.linalg as spla
12from fem_shell import ShellFEM, element_basis
14ch = json.load(open(sys.argv[1]))
15a = float(sys.argv[2]) if len(sys.argv) > 2 else 0.008
16k = int(sys.argv[3]) if len(sys.argv) > 3 else 14
17fem = ShellFEM(ch, a, log=lambda s: print(s))
18p = fem.p
19C, G, nu = fem.n_cells, fem.G, p["nu"]
20B = fem.B.reshape(G, 9, 48)
21# L: 9 derivative quantities -> [exx, eyy, gxy, kxx, kyy, 2kxy] at the reference state
22L = np.zeros((C, G, 6, 9))
23L[:, :, 0, 0] = fem.Zx; L[:, :, 0, 5] = 1.0
24L[:, :, 1, 1] = fem.Zy; L[:, :, 1, 8] = 1.0
25L[:, :, 2, 0] = fem.Zy; L[:, :, 2, 1] = fem.Zx; L[:, :, 2, 6] = 1.0; L[:, :, 2, 7] = 1.0
26L[:, :, 3, 2] = 1.0; L[:, :, 4, 3] = 1.0; L[:, :, 5, 4] = 2.0
27Q = np.array([[1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]])
28Dm = np.zeros((C, G, 6, 6))
29Dm[:, :, :3, :3] = fem.Cm[:, :, None, None] * Q
30Dm[:, :, 3:, 3:] = fem.Db[:, :, None, None] * Q
31LB = np.einsum("cgij,gjk->cgik", L, B) # (C, G, 6, 48)
32Ke = np.einsum("g,cgik,cgij,cgjl->ckl", fem.wq, LB, Dm, LB, optimize=True)
33rows = np.repeat(fem.cell_dofs, 48, axis=1).ravel()
34cols = np.tile(fem.cell_dofs, (1, 48)).ravel()
35K = sp.coo_matrix((Ke.ravel(), (rows, cols)), shape=(fem.n_dof, fem.n_dof)).tocsr()
36K = 0.5 * (K + K.T)
37# consistent mass
38N = element_basis(fem.sx, fem.sy, a)["N"] # (G, 16)
39Me16 = np.einsum("cg,gi,gj->cij", p["rho"] * fem.hq * fem.wq[None, :], N, N) # (C, 16, 16)
40Me = np.zeros((C, 48, 48))
41for b in range(3):
42 Me[:, 16 * b:16 * b + 16, 16 * b:16 * b + 16] = Me16
43Mc = sp.coo_matrix((Me.ravel(), (rows, cols)), shape=(fem.n_dof, fem.n_dof)).tocsr()
44Ml = sp.diags(fem.Mdiag)
45if fem.clamped:
46 fr = fem.free_dofs
47 K, Mc, Ml = K[fr][:, fr].tocsr(), Mc[fr][:, fr].tocsr(), sp.diags(fem.Mdiag[fr])
48 print(f"clamped at the mount: {len(fem.fixed_dofs)} unknowns fixed")
49print(f"total mass: consistent {Mc.sum():.5f}, lumped {fem.Mdiag[0:4 * fem.n_nodes:4].sum():.5f} (w block)")
50for name, M in [("lumped", Ml), ("consistent", Mc)]:
51 vals = spla.eigsh(K, k=k, M=M.tocsc(), sigma=-1e3, which="LM", return_eigenvectors=False)
52 f = np.sqrt(np.maximum(np.sort(vals), 0)) / (2 * np.pi)
53 print(f"{name:10s}: " + " ".join(f"{x:7.2f}" for x in f))
54# highest frequency with the lumped mass (stability)
55vals = spla.eigsh(K, k=1, M=Ml.tocsc(), which="LM", return_eigenvectors=False)
56print(f"omega_max (lumped) = {np.sqrt(vals[0]):.4e} rad/s")
57# which block sets the explicit step: bending-only and in-plane-only top frequencies
58nw = 4 * fem.n_nodes
59Kw = K[:nw, :nw]; Kp = K[nw:, nw:]
60for name, Kb, Mb in [("bending block", Kw, sp.diags(fem.Mdiag[:nw])), ("in-plane block", Kp, sp.diags(fem.Mdiag[nw:]))]:
61 v = spla.eigsh(Kb, k=1, M=Mb.tocsc(), which="LM", return_eigenvectors=False)
62 print(f"omega_max {name}: {np.sqrt(v[0]):.4e} rad/s")
63# and the bending block with the derivative-dof masses scaled up by 4 (rotational mass scaling)
64Ms = fem.Mdiag[:nw].copy(); Ms[1::4] *= 4; Ms[2::4] *= 4; Ms[3::4] *= 4
65v = spla.eigsh(Kw, k=1, M=sp.diags(Ms).tocsc(), which="LM", return_eigenvectors=False)
66print(f"omega_max bending with derivative masses x4: {np.sqrt(v[0]):.4e} rad/s")
67# in-plane block with its derivative-dof masses scaled, and the effect on the low modes of the whole system
68for sc in (4, 16, 64):
69 Ms = fem.Mdiag.copy()
70 for b in (1, 2):
71 base = nw * b
72 Ms[base + 1::4][: fem.n_nodes] *= sc; Ms[base + 2::4][: fem.n_nodes] *= sc; Ms[base + 3::4][: fem.n_nodes] *= sc
73 Msp = sp.diags(Ms)
74 v = spla.eigsh(Kp, k=1, M=sp.diags(Ms[nw:]).tocsc(), which="LM", return_eigenvectors=False)
75 low = spla.eigsh(K, k=14, M=Msp.tocsc(), sigma=-1e3, which="LM", return_eigenvectors=False)
76 fl = np.sqrt(np.maximum(np.sort(low), 0)) / (2 * np.pi)
77 print(f"in-plane derivative masses x{sc}: omega_max in-plane {np.sqrt(v[0]):.4e} rad/s; low modes " + " ".join(f"{x:6.2f}" for x in fl[6:]))
78# lowest in-plane modes with and without the derivative-mass scaling (in-plane block alone; 6 -> 3 rigid modes there)
79for sc in (1, 16):
80 Ms = fem.Mdiag[nw:].copy()
81 for b in (0, 1):
82 base = nw * b
83 Ms[base + 1::4][: fem.n_nodes] *= sc; Ms[base + 2::4][: fem.n_nodes] *= sc; Ms[base + 3::4][: fem.n_nodes] *= sc
84 low = spla.eigsh(Kp, k=9, M=sp.diags(Ms).tocsc(), sigma=-1e3, which="LM", return_eigenvectors=False)
85 fl = np.sqrt(np.maximum(np.sort(low), 0)) / (2 * np.pi)
86 print(f"in-plane modes, derivative masses x{sc}: " + " ".join(f"{x:7.1f}" for x in fl[3:]))