1// Stress-free initial shape of the sheet. The shape is a list of features
2// (bows, a twist, dents, creases, a random field) summed at the nodes, with
3// the mean plane removed. Only the slopes of the shape enter the physics, so
4// the mean is immaterial; the tilt is not, because the von Karman strains are
5// not invariant under rotation, and a tilted stress-free shape would add a
6// spurious membrane coupling. Every shape is therefore made plane-free here.
8import { type Grid, type PhysParams, insideNode, nodeWeight, nodeX, nodeY } from "./params.ts";
10export type FeatureType = "bow" | "twist" | "dent" | "crease" | "noise";
12/** Parabolic bow across the width (axis x) or the height (axis y); the offset in -1..1 moves the crest toward one edge. */
13export interface BowFeature { type: "bow"; axis: "x" | "y"; sag: number; offset: number }
14/** Saddle x*y with the given height at the corners (alternating in sign). */
15export interface TwistFeature { type: "twist"; height: number }
16/** Gaussian bump at (x, y) given as fractions of the width and the height, with the 1/e radius and the height. */
17export interface DentFeature { type: "dent"; x: number; y: number; radius: number; height: number }
18/** Fold along the line through (x, y) at `angle` degrees from the x axis, kinked by `fold` degrees (positive: ridge), rounded over `width`. */
19export interface CreaseFeature { type: "crease"; x: number; y: number; angle: number; fold: number; width: number }
20/** Seeded random field with Gaussian correlation of length `length` and the given RMS over the sheet. */
21export interface NoiseFeature { type: "noise"; seed: number; rms: number; length: number }
22export type ShapeFeature = BowFeature | TwistFeature | DentFeature | CreaseFeature | NoiseFeature;
24export interface FeatureField {
25 key: string;
26 label: string;
27 unit: string;
28 scale: number; // form value = SI value / scale
29 step: number;
30 options?: [string, string][]; // a choice rather than a number
31}
33export interface FeatureSpec {
34 label: string;
35 description: string;
36 fields: FeatureField[];
37}
39/** Presentation and encoding order of the fields of each feature type. */
40export const FEATURE_SPECS: Record<FeatureType, FeatureSpec> = {
41 bow: {
42 label: "Bow",
43 description: "Parabolic bow with zero sag along two opposite edges. The apex offset (-1 to 1) moves the crest toward one edge.",
44 fields: [
45 { key: "axis", label: "Across", unit: "", scale: 1, step: 1, options: [["x", "the width"], ["y", "the height"]] },
46 { key: "sag", label: "Sag", unit: "mm", scale: 1e-3, step: 0.5 },
47 { key: "offset", label: "Apex offset", unit: "", scale: 1, step: 0.1 },
48 ],
49 },
50 twist: {
51 label: "Twist",
52 description: "Saddle shape: two opposite corners up by the height, the other two down.",
53 fields: [{ key: "height", label: "Height at the corners", unit: "mm", scale: 1e-3, step: 0.5 }],
54 },
55 dent: {
56 label: "Dent",
57 description: "Gaussian bump; the radius is where the height has fallen to 1/e. A negative height is a dent toward the back.",
58 fields: [
59 { key: "x", label: "x (fraction of width)", unit: "", scale: 1, step: 0.05 },
60 { key: "y", label: "y (fraction of height)", unit: "", scale: 1, step: 0.05 },
61 { key: "radius", label: "Radius", unit: "mm", scale: 1e-3, step: 5 },
62 { key: "height", label: "Height", unit: "mm", scale: 1e-3, step: 0.5 },
63 ],
64 },
65 crease: {
66 label: "Crease",
67 description: "Fold along a line: the two sides tilt by half the fold angle each. A positive fold is a ridge, a negative one a valley; the width rounds the kink.",
68 fields: [
69 { key: "x", label: "x (fraction of width)", unit: "", scale: 1, step: 0.05 },
70 { key: "y", label: "y (fraction of height)", unit: "", scale: 1, step: 0.05 },
71 { key: "angle", label: "Line angle from x axis", unit: "deg", scale: 1, step: 15 },
72 { key: "fold", label: "Fold angle", unit: "deg", scale: 1, step: 0.5 },
73 { key: "width", label: "Rounding width", unit: "mm", scale: 1e-3, step: 5 },
74 ],
75 },
76 noise: {
77 label: "Random field",
78 description: "Smooth random undulation with Gaussian correlation; the same seed gives the same field on any grid.",
79 fields: [
80 { key: "seed", label: "Seed", unit: "", scale: 1, step: 1 },
81 { key: "rms", label: "RMS amplitude", unit: "mm", scale: 1e-3, step: 0.1 },
82 { key: "length", label: "Correlation length", unit: "mm", scale: 1e-3, step: 10 },
83 ],
84 },
85};
87export const FEATURE_TYPES = Object.keys(FEATURE_SPECS) as FeatureType[];
89/** A feature of the given type with sensible default values. */
90export function newFeature(type: FeatureType): ShapeFeature {
91 switch (type) {
92 case "bow": return { type, axis: "x", sag: 3e-3, offset: 0 };
93 case "twist": return { type, height: 1e-3 };
94 case "dent": return { type, x: 0.5, y: 0.5, radius: 0.03, height: 1e-3 };
95 case "crease": return { type, x: 0.5, y: 0.5, angle: 30, fold: 1, width: 0.015 };
96 case "noise": return { type, seed: 1, rms: 0.5e-3, length: 0.06 };
97 }
98}
100export function cloneFeatures(features: ShapeFeature[]): ShapeFeature[] {
101 return features.map((f) => ({ ...f }));
102}
104/** The shapes of the earlier interface: a parabolic bow across the width, the height, or both (dish). */
105export function legacyFeatures(bowShape: string, amplitude: number): ShapeFeature[] {
106 if (bowShape === "none" || amplitude === 0) return [];
107 if (bowShape === "height") return [{ type: "bow", axis: "y", sag: amplitude, offset: 0 }];
108 if (bowShape === "dish") return [{ type: "bow", axis: "x", sag: amplitude / 2, offset: 0 }, { type: "bow", axis: "y", sag: amplitude / 2, offset: 0 }];
109 return [{ type: "bow", axis: "x", sag: amplitude, offset: 0 }];
110}
112// ---------- Text encoding (URL) ----------
114function fmt(v: number): string {
115 return v === 0 ? "0" : String(+v.toPrecision(8));
116}
118/** Compact text form: `type:v1,v2,...;type:...` with the values in SI, or `none` for a flat sheet. */
119export function encodeFeatures(features: ShapeFeature[]): string {
120 if (features.length === 0) return "none";
121 return features
122 .map((f) => `${f.type}:${FEATURE_SPECS[f.type].fields.map((fld) => {
123 const v = (f as unknown as Record<string, number | string>)[fld.key];
124 return typeof v === "number" ? fmt(v) : String(v);
125 }).join(",")}`)
126 .join(";");
127}
129/** Inverse of encodeFeatures. Returns null when nothing in the text can be read; entries that cannot be read are skipped. */
130export function parseFeatures(text: string): ShapeFeature[] | null {
131 const s = text.trim();
132 if (s === "" || s === "none" || s === "flat") return [];
133 const out: ShapeFeature[] = [];
134 let any = false;
135 for (const item of s.split(";")) {
136 const colon = item.indexOf(":");
137 if (colon < 0) continue;
138 const type = item.slice(0, colon).trim() as FeatureType;
139 const spec = FEATURE_SPECS[type];
140 if (!spec) continue;
141 const vals = item.slice(colon + 1).split(",").map((v) => v.trim());
142 if (vals.length !== spec.fields.length) continue;
143 const f = newFeature(type) as unknown as Record<string, number | string>;
144 let ok = true;
145 spec.fields.forEach((fld, i) => {
146 if (fld.options) {
147 if (!fld.options.some(([v]) => v === vals[i])) ok = false;
148 else f[fld.key] = vals[i];
149 } else {
150 const num = parseFloat(vals[i]);
151 if (!isFinite(num)) ok = false;
152 else f[fld.key] = num;
153 }
154 });
155 if (!ok) continue;
156 out.push(f as unknown as ShapeFeature);
157 any = true;
158 }
159 return any ? out : null;
160}
162// ---------- Evaluation ----------
164/** Deterministic uniform generator on [0, 1) (mulberry32). */
165function mulberry32(seed: number): () => number {
166 let a = Math.floor(seed) >>> 0;
167 return () => {
168 a = (a + 0x6d2b79f5) >>> 0;
169 let t = a;
170 t = Math.imul(t ^ (t >>> 15), t | 1);
171 t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
172 return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
173 };
174}
176const NOISE_TERMS = 64;
178/**
179 * Adds the field of one feature to f (float64 per node, in meters). Nodes
180 * outside the sheet are filled too and zeroed at the end.
181 */
182function addFeature(f: Float64Array, feat: ShapeFeature, p: PhysParams, g: Grid): void {
183 const { nx, ny } = g;
184 if (feat.type === "noise") {
185 // Sum of random plane waves whose wavevectors are Gaussian with standard
186 // deviation sqrt(2)/L per component, which gives the covariance exp(-r^2/L^2).
187 const rnd = mulberry32(feat.seed);
188 const gauss = () => Math.sqrt(-2 * Math.log(1 - rnd())) * Math.cos(2 * Math.PI * rnd());
189 const L = Math.max(feat.length, 1e-4);
190 const terms = Array.from({ length: NOISE_TERMS }, () => ({ kx: (gauss() * Math.SQRT2) / L, ky: (gauss() * Math.SQRT2) / L, ph: 2 * Math.PI * rnd() }));
191 const raw = new Float64Array(g.n);
192 let sw = 0, s1 = 0, s2 = 0;
193 for (let j = 0; j < ny; j++) {
194 for (let i = 0; i < nx; i++) {
195 const x = nodeX(g, i), y = nodeY(g, j);
196 let v = 0;
197 for (const t of terms) v += Math.cos(t.kx * x + t.ky * y + t.ph);
198 raw[j * nx + i] = v;
199 const wgt = nodeWeight(g, i, j);
200 sw += wgt; s1 += wgt * v; s2 += wgt * v * v;
201 }
202 }
203 const mean = s1 / sw, rms = Math.sqrt(Math.max(0, s2 / sw - mean * mean));
204 const scale = rms > 0 ? feat.rms / rms : 0;
205 for (let k = 0; k < g.n; k++) f[k] += scale * (raw[k] - mean);
206 return;
207 }
208 for (let j = 0; j < ny; j++) {
209 for (let i = 0; i < nx; i++) {
210 const x = nodeX(g, i), y = nodeY(g, j);
211 const X = (2 * x) / p.Lx, Y = (2 * y) / p.Ly; // -1..1 over the bounding box
212 let v = 0;
213 switch (feat.type) {
214 case "bow": {
215 const s = feat.axis === "x" ? X : Y;
216 v = feat.sag * (1 - s * s) * (1 + feat.offset * s);
217 break;
218 }
219 case "twist":
220 v = feat.height * X * Y;
221 break;
222 case "dent": {
223 const dx = x - (feat.x - 0.5) * p.Lx, dy = y - (feat.y - 0.5) * p.Ly;
224 const r = Math.max(feat.radius, 1e-4);
225 v = feat.height * Math.exp(-(dx * dx + dy * dy) / (r * r));
226 break;
227 }
228 case "crease": {
229 const th = (feat.angle * Math.PI) / 180;
230 const d = -(x - (feat.x - 0.5) * p.Lx) * Math.sin(th) + (y - (feat.y - 0.5) * p.Ly) * Math.cos(th);
231 const s = Math.max(feat.width, 0);
232 v = -Math.tan((feat.fold * Math.PI) / 360) * (Math.sqrt(d * d + s * s) - s);
233 break;
234 }
235 }
236 f[j * nx + i] += v;
237 }
238 }
239}
241/** Subtracts the mass-weighted least-squares plane a + b x + c y from f over the sheet nodes. */
242export function removePlane(f: Float64Array, g: Grid): void {
243 let sw = 0, sx = 0, sy = 0, sxx = 0, sxy = 0, syy = 0, sf = 0, sxf = 0, syf = 0;
244 for (let j = 0; j < g.ny; j++) {
245 for (let i = 0; i < g.nx; i++) {
246 const wgt = nodeWeight(g, i, j);
247 if (wgt === 0) continue;
248 const x = nodeX(g, i), y = nodeY(g, j), v = f[j * g.nx + i];
249 sw += wgt; sx += wgt * x; sy += wgt * y;
250 sxx += wgt * x * x; sxy += wgt * x * y; syy += wgt * y * y;
251 sf += wgt * v; sxf += wgt * x * v; syf += wgt * y * v;
252 }
253 }
254 if (sw === 0) return;
255 // Normal equations of the 3-parameter fit, solved by Cramer's rule.
256 const det = sw * (sxx * syy - sxy * sxy) - sx * (sx * syy - sxy * sy) + sy * (sx * sxy - sxx * sy);
257 if (Math.abs(det) < 1e-300) return;
258 const a = (sf * (sxx * syy - sxy * sxy) - sx * (sxf * syy - sxy * syf) + sy * (sxf * sxy - sxx * syf)) / det;
259 const b = (sw * (sxf * syy - sxy * syf) - sf * (sx * syy - sxy * sy) + sy * (sx * syf - sxf * sy)) / det;
260 const c = (sw * (sxx * syf - sxf * sxy) - sx * (sx * syf - sxf * sy) + sf * (sx * sxy - sxx * sy)) / det;
261 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) f[j * g.nx + i] -= a + b * nodeX(g, i) + c * nodeY(g, j);
262}
264/** Stress-free initial shape w0 at every node: the sum of the features with the mean plane removed; zero outside the sheet. */
265export function initialShape(p: PhysParams, g: Grid): Float32Array<ArrayBuffer> {
266 const w0 = new Float32Array(g.n);
267 if (p.features.length === 0) return w0;
268 const f = new Float64Array(g.n);
269 for (const feat of p.features) addFeature(f, feat, p, g);
270 removePlane(f, g);
271 for (let j = 0; j < g.ny; j++) for (let i = 0; i < g.nx; i++) if (insideNode(g, i, j)) w0[j * g.nx + i] = f[j * g.nx + i];
272 return w0;
273}
275/** Peak-to-peak height and the largest slope of a shape over the sheet nodes. */
276export function shapeStats(w0: ArrayLike<number>, g: Grid): { peakToPeak: number; maxSlope: number } {
277 let lo = Infinity, hi = -Infinity, maxSlope = 0;
278 const at = (i: number, j: number) => w0[j * g.nx + i];
279 const slope = (i: number, j: number, di: number, dj: number, h: number) => {
280 const hasP = insideNode(g, i + di, j + dj), hasM = insideNode(g, i - di, j - dj);
281 if (hasP && hasM) return (at(i + di, j + dj) - at(i - di, j - dj)) / (2 * h);
282 if (hasP) return (at(i + di, j + dj) - at(i, j)) / h;
283 if (hasM) return (at(i, j) - at(i - di, j - dj)) / h;
284 return 0;
285 };
286 for (let j = 0; j < g.ny; j++) {
287 for (let i = 0; i < g.nx; i++) {
288 if (!insideNode(g, i, j)) continue;
289 const v = at(i, j);
290 lo = Math.min(lo, v); hi = Math.max(hi, v);
291 maxSlope = Math.max(maxSlope, Math.hypot(slope(i, j, 1, 0, g.hx), slope(i, j, 0, 1, g.hy)));
292 }
293 }
294 return { peakToPeak: isFinite(hi - lo) ? hi - lo : 0, maxSlope };
295}
297// ---------- Random configurations ----------
299/**
300 * A random feature list whose heights are of the order of `scale` (m): always
301 * a bow across one axis, sometimes one across the other, a twist, up to two
302 * dents, a crease, and a random field. `rnd` is a uniform generator on [0, 1).
303 */
304export function randomFeatures(p: PhysParams, scale: number, rnd: () => number = Math.random): ShapeFeature[] {
305 const U = (a: number, b: number) => a + (b - a) * rnd();
306 const sign = () => (rnd() < 0.5 ? -1 : 1);
307 const Lmin = Math.min(p.Lx, p.Ly);
308 const out: ShapeFeature[] = [];
309 const first: "x" | "y" = rnd() < 0.5 ? "x" : "y";
310 out.push({ type: "bow", axis: first, sag: sign() * U(0.3, 1) * scale, offset: U(-0.6, 0.6) });
311 if (rnd() < 0.5) out.push({ type: "bow", axis: first === "x" ? "y" : "x", sag: sign() * U(0.1, 0.5) * scale, offset: U(-0.6, 0.6) });
312 if (rnd() < 0.5) out.push({ type: "twist", height: sign() * U(0.1, 0.5) * scale });
313 const dents = rnd() < 0.3 ? 0 : rnd() < 0.6 ? 1 : 2;
314 for (let k = 0; k < dents; k++) out.push({ type: "dent", x: U(0.15, 0.85), y: U(0.15, 0.85), radius: U(0.1, 0.2) * Lmin, height: sign() * U(0.2, 0.7) * scale });
315 if (rnd() < 0.3) {
316 // A rise of a fraction of the scale over a quarter of the smaller side.
317 const slope = (U(0.2, 0.6) * scale) / (Lmin / 4);
318 out.push({ type: "crease", x: U(0.3, 0.7), y: U(0.3, 0.7), angle: U(0, 180), fold: sign() * ((2 * Math.atan(slope) * 180) / Math.PI), width: 0.05 * Lmin });
319 }
320 if (rnd() < 0.7) out.push({ type: "noise", seed: 1 + Math.floor(rnd() * 99999), rms: U(0.05, 0.25) * scale, length: U(0.1, 0.35) * Lmin });
321 // Round the values to what the form shows, so the list reads the same in the form and in a link.
322 for (const f of out) {
323 const rec = f as unknown as Record<string, number | string>;
324 for (const fld of FEATURE_SPECS[f.type].fields) {
325 if (fld.options || fld.key === "seed") continue;
326 const v = rec[fld.key] as number;
327 rec[fld.key] = +((v / fld.scale).toPrecision(4)) * fld.scale;
328 }
329 }
330 return out;
331}