Profile with hammer dents; full basis when the surface is not symmetric
2 changed files+58−20
solvers/vonkarman-modal-disk/__pycache__/vk_modal_disk.cpython-313.pycmodified+0−0View file
Binary file not shown.
solvers/vonkarman-modal-disk/vk_modal_disk.pymodified+58−20View file
@@ -28,7 +28,8 @@ modes needed for a given frequency reach (the coupling cost scales with the cube
2828 uses the unrestricted bases instead.
2929
3030 If the challenge has a shell-of-revolution geometry (cymbal-bow-v1: a paraboloidal bow z = Z(r);
31-cymbal-bell-v1: the bow plus a cosine-dome bell), the same solver applies the Marguerre shallow-shell equations, which add the initial shape w0 = Z to the von Karman
31+cymbal-bell-v1: the bow plus a cosine-dome bell; cymbal-hammered-v1: with raised-cosine dents, which break the
32+symmetry in y so the full modal basis is used), the same solver applies the Marguerre shallow-shell equations, which add the initial shape w0 = Z to the von Karman
3233 system: L(w + w0, Phi) on the right of the bending equation and 2 L(w0, w) in the compatibility equation.
3334 With G^p_j = int Psi_p L(w0, Phi_j) dA this gives a linear membrane stiffening 2 sum_p |c_p| G^p (G^p)^T,
3435 which is diagonalized once in the flat-mode basis to give the shell modes (exact in the Gautschi step,
@@ -65,9 +66,15 @@ class Profile:
6566 bell of cymbal-bell-v1. Uniform thickness along the normal. The hole of the bell challenge is not represented
6667 (the solver closes it; see the module docstring)."""
6768
68- def __init__(self, R: float, H: float, h: float, bell_radius: float = 0.0, bell_height: float = 0.0, hole_radius: float = 0.0):
69+ def __init__(self, R: float, H: float, h: float, bell_radius: float = 0.0, bell_height: float = 0.0, hole_radius: float = 0.0, dents=()):
6970 self.R, self.H, self.h = R, H, h
7071 self.rb, self.hb, self.rh = bell_radius, bell_height, hole_radius
72+ self.dents = [(float(d["x"]), float(d["y"]), float(d["radius"]), float(d["depth"])) for d in dents]
73+
74+ @property
75+ def symmetric_y(self) -> bool:
76+ """True when the surface is even in y (no dents): the axisymmetric bow and bell."""
77+ return len(self.dents) == 0
7178
7279 @staticmethod
7380 def from_challenge(ch: dict):
@@ -75,7 +82,7 @@ class Profile:
7582 if not g or g.get("type") != "shell-of-revolution":
7683 return None
7784 p = ch["parameters"]
78- return Profile(p["disk_radius"], p["bow_height"], p["h"], p.get("bell_radius", 0.0), p.get("bell_height", 0.0), p.get("hole_radius", 0.0))
85+ return Profile(p["disk_radius"], p["bow_height"], p["h"], p.get("bell_radius", 0.0), p.get("bell_height", 0.0), p.get("hole_radius", 0.0), g.get("dents", ()))
7986
8087 # ---- radial profile and its derivatives
8188 def _bell(self, r):
@@ -98,25 +105,51 @@ class Profile:
98105 zb, dzb, d2zb = self._bell(r)
99106 return z + zb, dz + dzb, d2z + d2zb
100107
108+ @staticmethod
109+ def _dent(rho, rd, depth):
110+ """A raised-cosine depression: value, first and second radial derivatives at distance rho from its center."""
111+ k = np.pi / rd
112+ inside = rho < rd
113+ z = np.where(inside, -0.5 * depth * (1.0 + np.cos(k * rho)), 0.0)
114+ dz = np.where(inside, 0.5 * depth * k * np.sin(k * rho), 0.0)
115+ d2z = np.where(inside, 0.5 * depth * k * k * np.cos(k * rho), 0.0)
116+ return z, dz, d2z
117+
118+ @staticmethod
119+ def _cartesian(x, y, dz, d2z):
120+ """Gradient and second derivatives of a radial function about the origin of (x, y), from Z'(r) and Z''(r)."""
121+ r = np.sqrt(x * x + y * y)
122+ rr = np.where(r > 0, r, 1.0)
123+ s = np.where(r > 0, dz / rr, d2z) # Z'/r -> Z'' on the axis
124+ cx, cy = x / rr, y / rr
125+ return s * x, s * y, d2z * cx * cx + s * cy * cy, d2z * cy * cy + s * cx * cx, (d2z - s) * cx * cy
126+
101127 def z(self, x, y):
102- return self.radial(np.sqrt(x * x + y * y))[0]
128+ x = np.asarray(x, dtype=float); y = np.asarray(y, dtype=float)
129+ z = self.radial(np.sqrt(x * x + y * y))[0]
130+ for xc, yc, rd, depth in self.dents:
131+ z = z + self._dent(np.hypot(x - xc, y - yc), rd, depth)[0]
132+ return z
103133
104134 def grad(self, x, y):
105- r = np.sqrt(x * x + y * y)
106- _, dz, _ = self.radial(r)
107- s = np.where(r > 0, dz / np.where(r > 0, r, 1.0), 0.0) # Z'/r, finite on the axis
108- return s * x, s * y
135+ x = np.asarray(x, dtype=float); y = np.asarray(y, dtype=float)
136+ _, dz, d2z = self.radial(np.sqrt(x * x + y * y))
137+ gx, gy, _, _, _ = self._cartesian(x, y, dz, d2z)
138+ for xc, yc, rd, depth in self.dents:
139+ _, ddz, dd2z = self._dent(np.hypot(x - xc, y - yc), rd, depth)
140+ ex, ey, _, _, _ = self._cartesian(x - xc, y - yc, ddz, dd2z)
141+ gx, gy = gx + ex, gy + ey
142+ return gx, gy
109143
110144 def curvatures(self, x, y):
111145 """Second derivatives (Z_xx, Z_yy, Z_xy) of the mid-surface height."""
112- r = np.sqrt(x * x + y * y)
113- _, dz, d2z = self.radial(r)
114- rr = np.where(r > 0, r, 1.0)
115- s = np.where(r > 0, dz / rr, d2z) # Z'/r -> Z'' on the axis
116- cx, cy = x / rr, y / rr
117- zxx = d2z * cx * cx + s * cy * cy
118- zyy = d2z * cy * cy + s * cx * cx
119- zxy = (d2z - s) * cx * cy
146+ x = np.asarray(x, dtype=float); y = np.asarray(y, dtype=float)
147+ _, dz, d2z = self.radial(np.sqrt(x * x + y * y))
148+ _, _, zxx, zyy, zxy = self._cartesian(x, y, dz, d2z)
149+ for xc, yc, rd, depth in self.dents:
150+ _, ddz, dd2z = self._dent(np.hypot(x - xc, y - yc), rd, depth)
151+ _, _, exx, eyy, exy = self._cartesian(x - xc, y - yc, ddz, dd2z)
152+ zxx, zyy, zxy = zxx + exx, zyy + eyy, zxy + exy
120153 return zxx, zyy, zxy
121154
122155 def normal(self, x, y):
@@ -130,8 +163,11 @@ class Profile:
130163 return x + 0.5 * self.h * nx, y + 0.5 * self.h * ny, self.z(x, y) + 0.5 * self.h * nz
131164
132165 def max_slope(self):
133- r = np.linspace(0, self.R, 20001)
134- return float(np.abs(self.radial(r)[1]).max())
166+ g = np.linspace(-self.R, self.R, 801)
167+ X, Y = np.meshgrid(g, g, indexing="ij")
168+ m = X * X + Y * Y <= self.R**2
169+ gx, gy = self.grad(X[m], Y[m])
170+ return float(np.sqrt(gx * gx + gy * gy).max())
135171
136172
137173 Paraboloid = Profile # former name
@@ -155,10 +191,12 @@ class VonKarmanDisk:
155191 ch2["parameters"] = p
156192 self.p = p
157193 self.log = log or (lambda s: None)
194+ self.profile = Profile.from_challenge(ch)
158195 if symmetric_y is None:
159- symmetric_y = p["strike_y"] == 0.0
196+ symmetric_y = p["strike_y"] == 0.0 and (self.profile is None or self.profile.symmetric_y)
197+ if symmetric_y and self.profile is not None and not self.profile.symmetric_y:
198+ raise ValueError("the dented surface is not symmetric in y; the symmetric basis cannot be used")
160199 self.symmetric_y = symmetric_y
161- self.profile = Profile.from_challenge(ch)
162200 t0 = time.perf_counter()
163201 self.plate = BFSDisk(ch2, spacing, n_w, symmetric_y=symmetric_y) # linear modes, modal forces, alpha, omega
164202 pl = self.plate