1"""The linearized plate as sparse matrices: the stiffness K (the Hessian of the
2discrete energy at the stress-free state) and the lumped mass M as a diagonal.
4A flat plate has one unknown per node (the transverse displacement w; the
5in-plane motion decouples and is neither struck nor heard). A plate with a
6profile couples w to the in-plane displacements u and v through the slopes of
7the profile, and has three unknowns per node, interleaved (u, v, w) so that
8the band stays narrow. The discretization is that of the full model
9(``full.kernel``): three-point second differences for the bending curvatures at
10the nodes, the four-point mixed difference at the cells, and bilinear cells
11with 2 x 2 Gauss integration for the membrane strains.
12"""
14from __future__ import annotations
16from dataclasses import dataclass
18import numpy as np
19import scipy.sparse as sp
21from ..plate import Plate
23G3 = 1 / np.sqrt(3.0)
26def shape_derivatives(gp: int, hx: float, hy: float) -> tuple[np.ndarray, np.ndarray]:
27 """Derivatives of the four bilinear shape functions at Gauss point ``gp`` (0..3).
28 The corners are ordered (ci, cj), (ci + 1, cj), (ci + 1, cj + 1), (ci, cj + 1)."""
29 xi = G3 if gp & 1 else -G3
30 eta = G3 if gp & 2 else -G3
31 dNdx = np.array([-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)]) / (2 * hx)
32 dNdy = np.array([-(1 - xi), -(1 + xi), 1 + xi, 1 - xi]) / (2 * hy)
33 return dNdx, dNdy
36def corners(a: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
37 """The values of a node array at the four corners of every cell, in the corner order above."""
38 return a[:-1, :-1], a[:-1, 1:], a[1:, 1:], a[1:, :-1]
42class LinearSystem:
43 K: sp.csr_matrix
44 """Symmetric stiffness, ``(n, n)``."""
45 mass: np.ndarray
46 """Diagonal mass of every unknown [kg], ``(n,)``."""
47 fields: int
48 """1 for a flat plate (w only), 3 for a plate with a profile (u, v, w interleaved)."""
49 dof: np.ndarray
50 """``(ny, nx)``: first unknown of each grid node, or -1 outside the plate."""
51 nodes: np.ndarray
52 """``(n_nodes,)`` flat grid index ``j * nx + i`` of each node, in unknown order."""
55 def n(self) -> int:
56 return len(self.mass)
59 def bandwidth(self) -> int:
60 """Half-bandwidth of K."""
61 coo = self.K.tocoo()
62 return int(np.abs(coo.row - coo.col).max()) if coo.nnz else 0
64 def w_dof(self) -> np.ndarray:
65 """``(ny, nx)``: the transverse unknown of each node, or -1 outside."""
66 return np.where(self.dof >= 0, self.dof + (self.fields - 1), -1)
68 def to_grid(self, x: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
69 """Scatters an unknown vector (or a stack ``(..., n)``) to grid arrays ``(u, v, w)``, zero outside the plate."""
70 x = np.asarray(x)
71 shape = x.shape[:-1] + (self.dof.shape[0] * self.dof.shape[1],)
72 u, v, w = np.zeros(shape), np.zeros(shape), np.zeros(shape)
73 if self.fields == 1:
74 w[..., self.nodes] = x
75 else:
76 u[..., self.nodes] = x[..., 0::3]
77 v[..., self.nodes] = x[..., 1::3]
78 w[..., self.nodes] = x[..., 2::3]
79 gs = x.shape[:-1] + self.dof.shape
80 return u.reshape(gs), v.reshape(gs), w.reshape(gs)
82 def from_grid(self, u: np.ndarray, v: np.ndarray, w: np.ndarray) -> np.ndarray:
83 """Gathers grid arrays into an unknown vector."""
84 if self.fields == 1:
85 return w.reshape(-1)[self.nodes].copy()
86 x = np.empty(self.n)
87 x[0::3] = u.reshape(-1)[self.nodes]
88 x[1::3] = v.reshape(-1)[self.nodes]
89 x[2::3] = w.reshape(-1)[self.nodes]
90 return x
93def assemble(pl: Plate, fields: int | None = None) -> LinearSystem:
94 """Assembles K and M for the plate. ``fields`` may force the three-field
95 system for a flat plate (the tests use it); by default a flat plate gets the
96 one-field system."""
97 g, m = pl.grid, pl.material
98 nx, ny, hx, hy = g.nx, g.ny, g.hx, g.hy
99 nf = fields if fields is not None else (1 if pl.is_flat else 3)
100 if nf not in (1, 3):
101 raise ValueError("fields must be 1 or 3")
102 nu = m.nu
104 # Number the plate nodes with the shorter grid axis innermost, so that the
105 # 13-point bending stencil, which reaches two rows away, gives the narrowest band.
106 inside = g.mask & (g.node_weight > 0)
107 if nx <= ny:
108 nodes = np.flatnonzero(inside.ravel())
109 else:
110 jj, ii = np.nonzero(inside.T) # (i, j) pairs in i-major order
111 nodes = ii * nx + jj
112 n_nodes = len(nodes)
113 n = n_nodes * nf
114 dof = -np.ones(g.n, dtype=np.int64)
115 dof[nodes] = np.arange(n_nodes) * nf
116 dof = dof.reshape(g.shape)
117 mass = np.repeat(pl.node_mass.reshape(-1)[nodes], nf)
118 wd = dof + (nf - 1) # w unknown of each node (garbage outside; masked below)
120 def stencil(rows_mask: np.ndarray, cols: list[np.ndarray], vals: list[float], n_rows: int) -> sp.csr_matrix:
121 """A sparse operator with one row per grid entity, nonzero where ``rows_mask``."""
122 r = np.flatnonzero(rows_mask.ravel())
123 data, ri, ci = [], [], []
124 for c, v in zip(cols, vals):
125 data.append(np.full(len(r), v))
126 ri.append(r)
127 ci.append(c.ravel()[r])
128 return sp.csr_matrix((np.concatenate(data), (np.concatenate(ri), np.concatenate(ci))), shape=(n_rows, n))
130 def shifted(a: np.ndarray, di: int, dj: int) -> np.ndarray:
131 """``a[j + dj, i + di]`` with -1 (an invalid column) outside the grid."""
132 out = -np.ones_like(a)
133 src = a[max(0, dj) : ny + min(0, dj), max(0, di) : nx + min(0, di)]
134 out[max(0, -dj) : ny + min(0, -dj), max(0, -di) : nx + min(0, -di)] = src
135 return out
137 ins = g.mask
138 left, right, down, up = g.neighbor_masks()
139 has_a = ins & left & right
140 has_b = ins & down & up
142 # ---- Bending: hx hy [ sum_k w_k D_k (a^2 + b^2 + 2 nu a b) + 2 (1 - nu) sum_c D_c c^2 ] ----
143 A = stencil(has_a, [shifted(wd, -1, 0), wd, shifted(wd, 1, 0)], [1 / hx**2, -2 / hx**2, 1 / hx**2], g.n)
144 B = stencil(has_b, [shifted(wd, 0, -1), wd, shifted(wd, 0, 1)], [1 / hy**2, -2 / hy**2, 1 / hy**2], g.n)
145 Fa = sp.diags((hx * hy * g.node_weight * m.bending_stiffness(pl.h)).ravel())
146 K = A.T @ Fa @ A + B.T @ Fa @ B + nu * (A.T @ Fa @ B + B.T @ Fa @ A)
147 cm = g.cell_mask
148 c00, c10, c11, c01 = corners(wd)
149 C = stencil(cm, [c00, c10, c01, c11], [1 / (hx * hy), -1 / (hx * hy), -1 / (hx * hy), 1 / (hx * hy)], cm.size)
150 Fc = sp.diags((2 * (1 - nu) * m.bending_stiffness(pl.cell_thickness) * hx * hy * cm).ravel())
151 K = K + C.T @ Fc @ C
153 # ---- Membrane, linearized about the profile (three-field system only) ----
154 if nf == 3:
155 q = hx * hy / 4
156 hc = pl.cell_thickness
157 Cm = sp.diags((q * m.membrane_stiffness(hc) * cm).ravel())
158 Gs = sp.diags((q * m.shear_stiffness(hc) * cm).ravel())
159 u_c = [c for c in corners(dof)]
160 v_c = [c + 1 for c in corners(dof)]
161 w_c = [c + 2 for c in corners(dof)]
162 w0_c = corners(pl.w0)
163 for gp in range(4):
164 dNdx, dNdy = shape_derivatives(gp, hx, hy)
165 w0x = sum(dNdx[b] * w0_c[b] for b in range(4))
166 w0y = sum(dNdy[b] * w0_c[b] for b in range(4))
167 # Linear strains: exx = u_x + w0_x w_x, eyy = v_y + w0_y w_y, gxy = u_y + v_x + w0_x w_y + w0_y w_x.
168 Bxx = var_stencil(cm, u_c + w_c, [np.full(cm.shape, dNdx[b]) for b in range(4)] + [w0x * dNdx[b] for b in range(4)], n)
169 Byy = var_stencil(cm, v_c + w_c, [np.full(cm.shape, dNdy[b]) for b in range(4)] + [w0y * dNdy[b] for b in range(4)], n)
170 Bxy = var_stencil(cm, u_c + v_c + w_c, [np.full(cm.shape, dNdy[b]) for b in range(4)] + [np.full(cm.shape, dNdx[b]) for b in range(4)] + [w0x * dNdy[b] + w0y * dNdx[b] for b in range(4)], n)
171 K = K + Bxx.T @ Cm @ Bxx + Byy.T @ Cm @ Byy + nu * (Bxx.T @ Cm @ Byy + Byy.T @ Cm @ Bxx) + Bxy.T @ Gs @ Bxy
173 K = sp.csr_matrix(K)
174 K = (K + K.T) * 0.5
175 K.eliminate_zeros()
176 K.sort_indices()
177 return LinearSystem(K, mass, nf, dof, nodes)
180def var_stencil(rows_mask: np.ndarray, cols: list[np.ndarray], vals: list[np.ndarray], n: int) -> sp.csr_matrix:
181 """A sparse operator with one row per cell and per-cell coefficients."""
182 r = np.flatnonzero(rows_mask.ravel())
183 data = np.concatenate([v.ravel()[r] for v in vals])
184 ri = np.concatenate([r] * len(cols))
185 ci = np.concatenate([c.ravel()[r] for c in cols])
186 return sp.csr_matrix((data, (ri, ci)), shape=(rows_mask.size, n))