"""The linearized plate as sparse matrices: the stiffness K (the Hessian of the discrete energy at the stress-free state) and the lumped mass M as a diagonal. A flat plate has one unknown per node (the transverse displacement w; the in-plane motion decouples and is neither struck nor heard). A plate with a profile couples w to the in-plane displacements u and v through the slopes of the profile, and has three unknowns per node, interleaved (u, v, w) so that the band stays narrow. The discretization is that of the full model (``full.kernel``): three-point second differences for the bending curvatures at the nodes, the four-point mixed difference at the cells, and bilinear cells with 2 x 2 Gauss integration for the membrane strains. """ from __future__ import annotations from dataclasses import dataclass import numpy as np import scipy.sparse as sp from ..plate import Plate G3 = 1 / np.sqrt(3.0) def shape_derivatives(gp: int, hx: float, hy: float) -> tuple[np.ndarray, np.ndarray]: """Derivatives of the four bilinear shape functions at Gauss point ``gp`` (0..3). The corners are ordered (ci, cj), (ci + 1, cj), (ci + 1, cj + 1), (ci, cj + 1).""" xi = G3 if gp & 1 else -G3 eta = G3 if gp & 2 else -G3 dNdx = np.array([-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)]) / (2 * hx) dNdy = np.array([-(1 - xi), -(1 + xi), 1 + xi, 1 - xi]) / (2 * hy) return dNdx, dNdy def corners(a: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """The values of a node array at the four corners of every cell, in the corner order above.""" return a[:-1, :-1], a[:-1, 1:], a[1:, 1:], a[1:, :-1] @dataclass class LinearSystem: K: sp.csr_matrix """Symmetric stiffness, ``(n, n)``.""" mass: np.ndarray """Diagonal mass of every unknown [kg], ``(n,)``.""" fields: int """1 for a flat plate (w only), 3 for a plate with a profile (u, v, w interleaved).""" dof: np.ndarray """``(ny, nx)``: first unknown of each grid node, or -1 outside the plate.""" nodes: np.ndarray """``(n_nodes,)`` flat grid index ``j * nx + i`` of each node, in unknown order.""" @property def n(self) -> int: return len(self.mass) @property def bandwidth(self) -> int: """Half-bandwidth of K.""" coo = self.K.tocoo() return int(np.abs(coo.row - coo.col).max()) if coo.nnz else 0 def w_dof(self) -> np.ndarray: """``(ny, nx)``: the transverse unknown of each node, or -1 outside.""" return np.where(self.dof >= 0, self.dof + (self.fields - 1), -1) def to_grid(self, x: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Scatters an unknown vector (or a stack ``(..., n)``) to grid arrays ``(u, v, w)``, zero outside the plate.""" x = np.asarray(x) shape = x.shape[:-1] + (self.dof.shape[0] * self.dof.shape[1],) u, v, w = np.zeros(shape), np.zeros(shape), np.zeros(shape) if self.fields == 1: w[..., self.nodes] = x else: u[..., self.nodes] = x[..., 0::3] v[..., self.nodes] = x[..., 1::3] w[..., self.nodes] = x[..., 2::3] gs = x.shape[:-1] + self.dof.shape return u.reshape(gs), v.reshape(gs), w.reshape(gs) def from_grid(self, u: np.ndarray, v: np.ndarray, w: np.ndarray) -> np.ndarray: """Gathers grid arrays into an unknown vector.""" if self.fields == 1: return w.reshape(-1)[self.nodes].copy() x = np.empty(self.n) x[0::3] = u.reshape(-1)[self.nodes] x[1::3] = v.reshape(-1)[self.nodes] x[2::3] = w.reshape(-1)[self.nodes] return x def assemble(pl: Plate, fields: int | None = None) -> LinearSystem: """Assembles K and M for the plate. ``fields`` may force the three-field system for a flat plate (the tests use it); by default a flat plate gets the one-field system.""" g, m = pl.grid, pl.material nx, ny, hx, hy = g.nx, g.ny, g.hx, g.hy nf = fields if fields is not None else (1 if pl.is_flat else 3) if nf not in (1, 3): raise ValueError("fields must be 1 or 3") nu = m.nu # Number the plate nodes with the shorter grid axis innermost, so that the # 13-point bending stencil, which reaches two rows away, gives the narrowest band. inside = g.mask & (g.node_weight > 0) if nx <= ny: nodes = np.flatnonzero(inside.ravel()) else: jj, ii = np.nonzero(inside.T) # (i, j) pairs in i-major order nodes = ii * nx + jj n_nodes = len(nodes) n = n_nodes * nf dof = -np.ones(g.n, dtype=np.int64) dof[nodes] = np.arange(n_nodes) * nf dof = dof.reshape(g.shape) mass = np.repeat(pl.node_mass.reshape(-1)[nodes], nf) wd = dof + (nf - 1) # w unknown of each node (garbage outside; masked below) def stencil(rows_mask: np.ndarray, cols: list[np.ndarray], vals: list[float], n_rows: int) -> sp.csr_matrix: """A sparse operator with one row per grid entity, nonzero where ``rows_mask``.""" r = np.flatnonzero(rows_mask.ravel()) data, ri, ci = [], [], [] for c, v in zip(cols, vals): data.append(np.full(len(r), v)) ri.append(r) ci.append(c.ravel()[r]) return sp.csr_matrix((np.concatenate(data), (np.concatenate(ri), np.concatenate(ci))), shape=(n_rows, n)) def shifted(a: np.ndarray, di: int, dj: int) -> np.ndarray: """``a[j + dj, i + di]`` with -1 (an invalid column) outside the grid.""" out = -np.ones_like(a) src = a[max(0, dj) : ny + min(0, dj), max(0, di) : nx + min(0, di)] out[max(0, -dj) : ny + min(0, -dj), max(0, -di) : nx + min(0, -di)] = src return out ins = g.mask left, right, down, up = g.neighbor_masks() has_a = ins & left & right has_b = ins & down & up # ---- 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 ] ---- A = stencil(has_a, [shifted(wd, -1, 0), wd, shifted(wd, 1, 0)], [1 / hx**2, -2 / hx**2, 1 / hx**2], g.n) B = stencil(has_b, [shifted(wd, 0, -1), wd, shifted(wd, 0, 1)], [1 / hy**2, -2 / hy**2, 1 / hy**2], g.n) Fa = sp.diags((hx * hy * g.node_weight * m.bending_stiffness(pl.h)).ravel()) K = A.T @ Fa @ A + B.T @ Fa @ B + nu * (A.T @ Fa @ B + B.T @ Fa @ A) cm = g.cell_mask c00, c10, c11, c01 = corners(wd) C = stencil(cm, [c00, c10, c01, c11], [1 / (hx * hy), -1 / (hx * hy), -1 / (hx * hy), 1 / (hx * hy)], cm.size) Fc = sp.diags((2 * (1 - nu) * m.bending_stiffness(pl.cell_thickness) * hx * hy * cm).ravel()) K = K + C.T @ Fc @ C # ---- Membrane, linearized about the profile (three-field system only) ---- if nf == 3: q = hx * hy / 4 hc = pl.cell_thickness Cm = sp.diags((q * m.membrane_stiffness(hc) * cm).ravel()) Gs = sp.diags((q * m.shear_stiffness(hc) * cm).ravel()) u_c = [c for c in corners(dof)] v_c = [c + 1 for c in corners(dof)] w_c = [c + 2 for c in corners(dof)] w0_c = corners(pl.w0) for gp in range(4): dNdx, dNdy = shape_derivatives(gp, hx, hy) w0x = sum(dNdx[b] * w0_c[b] for b in range(4)) w0y = sum(dNdy[b] * w0_c[b] for b in range(4)) # 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. 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) 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) 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) K = K + Bxx.T @ Cm @ Bxx + Byy.T @ Cm @ Byy + nu * (Bxx.T @ Cm @ Byy + Byy.T @ Cm @ Bxx) + Bxy.T @ Gs @ Bxy K = sp.csr_matrix(K) K = (K + K.T) * 0.5 K.eliminate_zeros() K.sort_indices() return LinearSystem(K, mass, nf, dof, nodes) def var_stencil(rows_mask: np.ndarray, cols: list[np.ndarray], vals: list[np.ndarray], n: int) -> sp.csr_matrix: """A sparse operator with one row per cell and per-cell coefficients.""" r = np.flatnonzero(rows_mask.ravel()) data = np.concatenate([v.ravel()[r] for v in vals]) ri = np.concatenate([r] * len(cols)) ci = np.concatenate([c.ravel()[r] for c in cols]) return sp.csr_matrix((data, (ri, ci)), shape=(rows_mask.size, n))