1"""The discrete potential energy of the plate in NumPy (float64), whose negative
2gradient is the force of the full model's kernel, and the largest frequency of
3the linearized system, which sets the time step. Neither is used in a time loop;
4they serve the tests and the set-up of a run."""
6from __future__ import annotations
8import numpy as np
9import scipy.linalg
11from ..grid import Grid
12from ..linear.assemble import assemble, corners, shape_derivatives
13from ..material import Material
14from ..plate import Plate
17def potential_energy(pl: Plate, u: np.ndarray, v: np.ndarray, w: np.ndarray) -> float:
18 """Bending plus membrane energy [J] of the displacement fields ``(u, v, w)``
19 ``(ny, nx)``, measured from the stress-free shape."""
20 g, m = pl.grid, pl.material
21 hx, hy, nu = g.hx, g.hy, m.nu
22 left, right, down, up = g.neighbor_masks()
23 wp = np.pad(w, 1)
24 a = np.where(left & right, (wp[1:-1, 2:] - 2 * w + wp[1:-1, :-2]) / hx**2, 0.0)
25 b = np.where(down & up, (wp[2:, 1:-1] - 2 * w + wp[:-2, 1:-1]) / hy**2, 0.0)
26 Ub = (g.node_weight * m.bending_stiffness(pl.h) * (a * a + b * b + 2 * nu * a * b)).sum()
27 cc = (w[1:, 1:] - w[1:, :-1] - w[:-1, 1:] + w[:-1, :-1]) / (hx * hy)
28 Ub += (2 * (1 - nu) * m.bending_stiffness(pl.cell_thickness) * cc * cc * g.cell_mask).sum()
29 Ub *= 0.5 * hx * hy
31 hc = pl.cell_thickness
32 C, Gs = m.membrane_stiffness(hc), m.shear_stiffness(hc)
33 Uc, Vc, Wc, W0c = corners(u), corners(v), corners(w), corners(pl.w0)
34 Um = 0.0
35 q = hx * hy / 4
36 for gp in range(4):
37 dNdx, dNdy = shape_derivatives(gp, hx, hy)
38 ux = sum(dNdx[k] * Uc[k] for k in range(4))
39 uy = sum(dNdy[k] * Uc[k] for k in range(4))
40 vx = sum(dNdx[k] * Vc[k] for k in range(4))
41 vy = sum(dNdy[k] * Vc[k] for k in range(4))
42 dwx = sum(dNdx[k] * Wc[k] for k in range(4))
43 dwy = sum(dNdy[k] * Wc[k] for k in range(4))
44 w0x = sum(dNdx[k] * W0c[k] for k in range(4))
45 w0y = sum(dNdy[k] * W0c[k] for k in range(4))
46 exx = ux + dwx * (w0x + 0.5 * dwx)
47 eyy = vy + dwy * (w0y + 0.5 * dwy)
48 gxy = uy + vx + w0x * dwy + w0y * dwx + dwx * dwy
49 Um += (q * 0.5 * (C * (exx * exx + eyy * eyy + 2 * nu * exx * eyy) + Gs * gxy * gxy) * g.cell_mask).sum()
50 return float(Ub + Um)
53def max_angular_frequency(g: Grid, m: Material, h: float, n_small: int = 12) -> float:
54 """Largest angular frequency of the linearized discrete system [rad/s], from
55 the three-field stiffness of a small uniform flat plate with the same
56 spacing and the thickness ``h`` (the largest of the plate, which bounds the
57 bending frequency; the in-plane wave speed does not depend on the
58 thickness). The extreme mode is a grid-scale mode, so neither the plate
59 size nor its outline matters."""
60 gs = Grid(n_small, n_small, g.hx, g.hy, np.ones((n_small, n_small), dtype=bool))
61 sys = assemble(Plate.uniform(gs, m, h), fields=3)
62 s = 1 / np.sqrt(sys.mass)
63 A = sys.K.multiply(s[:, None]).multiply(s[None, :]).toarray()
64 lam = scipy.linalg.eigvalsh(A, subset_by_index=(sys.n - 1, sys.n - 1))[0]
65 return float(np.sqrt(max(lam, 0.0)))