Clamped mount: fixed unknowns at the inner staircase edge in the FEM, free-unknown eigenproblem in the modal solver
3 changed files+37−11
solvers/hermite-fem/check_modes.pymodified+4−0View file
@@ -42,6 +42,10 @@ for b in range(3):
4242 Me[:, 16 * b:16 * b + 16, 16 * b:16 * b + 16] = Me16
4343 Mc = sp.coo_matrix((Me.ravel(), (rows, cols)), shape=(fem.n_dof, fem.n_dof)).tocsr()
4444 Ml = sp.diags(fem.Mdiag)
45+if 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")
4549 print(f"total mass: consistent {Mc.sum():.5f}, lumped {fem.Mdiag[0:4 * fem.n_nodes:4].sum():.5f} (w block)")
4650 for name, M in [("lumped", Ml), ("consistent", Mc)]:
4751 vals = spla.eigsh(K, k=k, M=M.tocsc(), sigma=-1e3, which="LM", return_eigenvectors=False)
solvers/hermite-fem/fem_shell.pymodified+16−4View file
@@ -43,7 +43,7 @@ class Profile:
4343
4444 def __init__(self, p: dict):
4545 self.R = p["disk_radius"]
46- self.rh = p["hole_radius"]
46+ self.rh = p.get("mount_radius", p.get("hole_radius", 0.0)) # inner radius of the shell (the clamped mount)
4747 self.H = p.get("bow_height", 0.0)
4848 self.rb = p.get("bell_radius", 0.0)
4949 self.hb = p.get("bell_height", 0.0)
@@ -137,6 +137,7 @@ class ShellFEM:
137137 keep = (rc < R) & (rc > rh)
138138 self.cells = np.stack([ii.ravel()[keep], jj.ravel()[keep]], axis=1) # lower-left node index (i, j)
139139 self.n_cells = len(self.cells)
140+ inner_cells = set(map(tuple, np.stack([ii.ravel()[rc <= rh], jj.ravel()[rc <= rh]], axis=1).tolist())) # cells inside the mount
140141 # node map
141142 node_ids = {}
142143 cell_nodes = np.zeros((self.n_cells, 4), dtype=np.int64)
@@ -152,6 +153,15 @@ class ShellFEM:
152153 self.node_xy[k] = (i * a, j * a)
153154 self.node_ids = node_ids
154155 self.n_dof = 12 * self.n_nodes # blocks w, u, v of 4 n_nodes each
156+ # the clamped inner edge: every node that touches a cell inside the mount has all its unknowns fixed at zero
157+ # (displacements and the slopes of w), which is the clamp of MODEL.md on the staircase inner boundary
158+ fixed_nodes = np.array([k for (i, j), k in node_ids.items() if any((i + di, j + dj) in inner_cells for di in (-1, 0) for dj in (-1, 0))], dtype=np.int64)
159+ self.fixed_dofs = np.concatenate([b * 4 * self.n_nodes + 4 * fixed_nodes[:, None] + np.arange(4)[None, :] for b in range(3)]).ravel() if len(fixed_nodes) else np.zeros(0, dtype=np.int64)
160+ self.free_dofs = np.setdiff1d(np.arange(self.n_dof), self.fixed_dofs)
161+ self.clamped = len(fixed_nodes) > 0
162+ if self.clamped and no_drift:
163+ no_drift = False # the clamp removes the rigid motions; no reaction or projection is needed
164+ self.no_drift = False
155165 loc = (4 * cell_nodes[:, :, None] + np.arange(4)[None, None, :]).reshape(self.n_cells, 16)
156166 self.cell_dofs = np.concatenate([loc, loc + 4 * self.n_nodes, loc + 8 * self.n_nodes], axis=1) # (C, 48): w | u | v
157167 self.cell_dofs_w = loc
@@ -263,8 +273,8 @@ class ShellFEM:
263273 d = np.sqrt((L[0] - gx) ** 2 + (L[1] - gy) ** 2 + (L[2] - gz) ** 2) # (C, G)
264274 self.ray_c = (p["rho_air"] / (2 * np.pi)) * self.wq[None, :] / d
265275 self.ray_tau = d / p["c_air"]
266- self.log(f" mesh: {self.n_cells} cells, {self.n_nodes} nodes, {self.n_dof} unknowns, spacing {a * 1e3:g} mm, mass {self.mass_total:.4f} kg, "
267- f"footprint integrals {[round(s['integral'], 4) for s in self.strikes]}")
276+ self.log(f" mesh: {self.n_cells} cells, {self.n_nodes} nodes, {self.n_dof} unknowns ({len(self.fixed_dofs)} fixed at the mount), spacing {a * 1e3:g} mm, "
277+ f"mass {self.mass_total:.4f} kg, footprint integrals {[round(s['integral'], 4) for s in self.strikes]}")
268278
269279 # ------------------------------------------------------------------------------------------------------------
270280 def assemble_linear(self):
@@ -345,7 +355,9 @@ class ShellFEM:
345355 f.index_add_(0, cell_dofs.view(-1), fe.view(-1))
346356 return f
347357
348- ops = dict(cell_dofs=cell_dofs, cell_dofs_w=T64i(self.cell_dofs_w), NwT=T64(self.Nw.T.copy()), invM=T64(1.0 / self.Mdiag), T64=T64, T64i=T64i)
358+ invM_np = 1.0 / self.Mdiag
359+ invM_np[self.fixed_dofs] = 0.0 # clamped unknowns never accelerate
360+ ops = dict(cell_dofs=cell_dofs, cell_dofs_w=T64i(self.cell_dofs_w), NwT=T64(self.Nw.T.copy()), invM=T64(invM_np), T64=T64, T64i=T64i)
349361 return force, ops
350362
351363 def field_interpolator(self, xs, ys):
solvers/hermite-fem/modal_shell.pymodified+17−7View file
@@ -47,7 +47,10 @@ def sliced_modes(K, Mdiag, rigid, f_max: float, k_slice: int = 250, log=print, p
4747 factorizing K - sigma M once (the low-rank term is handled with the Woodbury identity)."""
4848 n = K.shape[0]
4949 M = sp.diags(Mdiag).tocsc()
50- U = (Mdiag[:, None] * rigid.T) * math.sqrt(penalty) # (n, 6): penalty term is U U^T
50+ if rigid is None: # a clamped shell has no rigid fields: plain shift-invert
51+ rigid = np.zeros((0, n))
52+ U = (Mdiag[:, None] * rigid.T) * math.sqrt(penalty) # (n, r): penalty term is U U^T
53+ nr = U.shape[1]
5154 Kop = spla.LinearOperator((n, n), matvec=lambda x: K @ x + U @ (U.T @ x), dtype=np.float64)
5255 target = (2 * math.pi * f_max) ** 2
5356 vals_all, vecs_all = [], []
@@ -58,8 +61,8 @@ def sliced_modes(K, Mdiag, rigid, f_max: float, k_slice: int = 250, log=print, p
5861 while True:
5962 t0 = time.perf_counter()
6063 lu = spla.splu((K - sigma * M).tocsc())
61- AiU = lu.solve(U) # (n, 6)
62- S = np.linalg.inv(np.eye(6) + U.T @ AiU)
64+ AiU = lu.solve(U) if nr else np.zeros((n, 0)) # (n, r)
65+ S = np.linalg.inv(np.eye(nr) + U.T @ AiU) if nr else np.zeros((0, 0))
6366
6467 def opinv(x, lu=lu, AiU=AiU, S=S):
6568 y = lu.solve(x)
@@ -90,8 +93,8 @@ def sliced_modes(K, Mdiag, rigid, f_max: float, k_slice: int = 250, log=print, p
9093 sel = (vals <= target) & (vals < 0.5 * penalty)
9194 vals, vecs = vals[sel], vecs[:, sel]
9295 # rigid content check and clean-up of the near-null modes of the flat geometry (exact zeros the penalty removed)
93- rc = np.abs(rigid @ (Mdiag[:, None] * vecs)).max()
94- log(f" {len(vals)} modes below {f_max:g} Hz; largest rigid component {rc:.1e}")
96+ rc = np.abs(rigid @ (Mdiag[:, None] * vecs)).max() if nr else 0.0
97+ log(f" {len(vals)} modes below {f_max:g} Hz" + (f"; largest rigid component {rc:.1e}" if nr else " (clamped: no rigid fields)"))
9598 return np.sqrt(np.maximum(vals, 0.0)), vecs
9699
97100
@@ -107,14 +110,21 @@ class ModalShell:
107110 cache_dir.mkdir(exist_ok=True)
108111 import hashlib
109112
110- geom_key = hashlib.sha1(json.dumps({k: v for k, v in fem.ch["parameters"].items() if k in ("disk_radius", "hole_radius", "bow_height", "bell_radius", "bell_height", "thickness", "thickness_center", "thickness_rim", "E", "nu", "rho")}, sort_keys=True).encode()).hexdigest()[:8]
113+ geom_key = hashlib.sha1(json.dumps({k: v for k, v in fem.ch["parameters"].items() if k in ("disk_radius", "hole_radius", "mount_radius", "bow_height", "bell_radius", "bell_height", "thickness", "thickness_center", "thickness_rim", "E", "nu", "rho")}, sort_keys=True).encode()).hexdigest()[:8]
111114 cache = cache_dir / f"modes-{fem.ch['challenge_id']}-g{fem.a * 1e3:g}-f{f_max / 1e3:g}k-s{fem.inplane_scale:g}-{geom_key}.npz"
112115 if cache.exists():
113116 z = np.load(cache)
114117 self.omega, self.Phi = z["omega"], z["Phi"]
115118 self.log(f" modes loaded from {cache.name}")
116119 else:
117- self.omega, self.Phi = sliced_modes(K, fem.Mdiag, fem.rigid, f_max, log=log)
120+ if fem.clamped: # eigenproblem on the free unknowns only; the modes are zero at the mount
121+ fr = fem.free_dofs
122+ om, Pf = sliced_modes(K[fr][:, fr].tocsr(), fem.Mdiag[fr], None, f_max, log=log)
123+ Phi = np.zeros((fem.n_dof, len(om)))
124+ Phi[fr] = Pf
125+ self.omega, self.Phi = om, Phi
126+ else:
127+ self.omega, self.Phi = sliced_modes(K, fem.Mdiag, fem.rigid, f_max, log=log)
118128 np.savez(cache, omega=self.omega, Phi=self.Phi)
119129 self.n_modes = len(self.omega)
120130 self.log(f" {self.n_modes} modes in {time.perf_counter() - t0:.1f}s: f = {np.round(self.omega[:5] / 2 / np.pi, 2)} ... {self.omega[-1] / 2 / np.pi:.0f} Hz")