1/**
2 * Inverse metric quantities V_theta, V_phi of a surface embedding X=(x,y,z)
3 * (evolving_surface/notes/algos.tex Algorithm 2 / SurfaceDiffOperator.
4 * _precompute_metric_quantities, clear_denominators=False branch): six grid
5 * scalar fields depending only on the geometry, used by the surface
6 * Laplace-Beltrami operator (Algorithm 3) to contract a field's theta/phi
7 * derivatives into a tangential gradient/divergence.
8 *
9 * g_tt = Xt.Xt, g_tp = Xt.Xp, g_pp = Xp.Xp (first fundamental form)
10 * det = g_tt*g_pp - g_tp^2
11 * V_theta = ( g_pp*Xt - g_tp*Xp ) / det
12 * V_phi = ( g_tt*Xp - g_tp*Xt ) / det
13 */
15export interface MetricFields {
16 /** V_theta, Cartesian components, npts each. */
17 Vtx: Float32Array;
18 Vty: Float32Array;
19 Vtz: Float32Array;
20 /** V_phi, Cartesian components, npts each. */
21 Vpx: Float32Array;
22 Vpy: Float32Array;
23 Vpz: Float32Array;
24}
26/**
27 * Xt/Xp (etc) are the theta/phi derivatives of each Cartesian embedding
28 * component, grid space, npts each -- the tangent vectors X_theta, X_phi of
29 * algos.tex Sec 4.1, one component per array.
30 */
31export function computeMetric(
32 npts: number,
33 Xt: Float32Array,
34 Xp: Float32Array,
35 Yt: Float32Array,
36 Yp: Float32Array,
37 Zt: Float32Array,
38 Zp: Float32Array,
39): MetricFields {
40 const Vtx = new Float32Array(npts);
41 const Vty = new Float32Array(npts);
42 const Vtz = new Float32Array(npts);
43 const Vpx = new Float32Array(npts);
44 const Vpy = new Float32Array(npts);
45 const Vpz = new Float32Array(npts);
47 for (let i = 0; i < npts; i++) {
48 const xt = Xt[i];
49 const xp = Xp[i];
50 const yt = Yt[i];
51 const yp = Yp[i];
52 const zt = Zt[i];
53 const zp = Zp[i];
55 const gtt = xt * xt + yt * yt + zt * zt;
56 const gtp = xt * xp + yt * yp + zt * zp;
57 const gpp = xp * xp + yp * yp + zp * zp;
58 const det = gtt * gpp - gtp * gtp;
60 Vtx[i] = (gpp * xt - gtp * xp) / det;
61 Vty[i] = (gpp * yt - gtp * yp) / det;
62 Vtz[i] = (gpp * zt - gtp * zp) / det;
63 Vpx[i] = (gtt * xp - gtp * xt) / det;
64 Vpy[i] = (gtt * yp - gtp * yt) / det;
65 Vpz[i] = (gtt * zp - gtp * zt) / det;
66 }
68 return { Vtx, Vty, Vtz, Vpx, Vpy, Vpz };
69}
71/**
72 * Flux-form metric weights p1, p2, q2, r of the six-transform Laplace-Beltrami
73 * scheme (docs/reduced-transforms.md Sec 3). Built from the
74 * *sin-weighted* theta tangent sin(theta)*X_theta — the undivided synthesis of
75 * the alpha shift, DerivPlan.sinDtheta — and X_phi, both smooth on the sphere:
76 *
77 * gtt~ = |sin(theta) X_theta|^2 (= sin^2(theta) g_tt)
78 * gtp~ = (sin(theta) X_theta).X_phi (= sin(theta) g_tp)
79 * gpp = |X_phi|^2
80 * D = sqrt(gtt~ gpp - gtp~^2) (= sin^2(theta) sqrt(det g) / sin(theta)
81 * = J sin^2(theta), with J = sqrt(det g)/sin(theta))
82 *
83 * p1 = gpp / D, p2 = -gtp~ / D, q2 = gtt~ / D, r = 1 / D.
84 *
85 * With these, for A = sin(theta) dtheta(u) and B = dphi(u), the two fluxes
86 *
87 * P = p1*A + p2*B, Qtilde = p2*A + q2*B
88 *
89 * equal sqrt(det g) g^{theta j} u_j and sin(theta) sqrt(det g) g^{phi j} u_j —
90 * both smooth on the sphere — and Delta_Gamma u = r * (sin(theta) dtheta(P) +
91 * dphi(Qtilde)). p1, p2, q2 are bounded (the sin^2 in D cancels against the
92 * vanishing numerators); r ~ 1/sin^2(theta) is finite at the Gauss nodes and
93 * is the scheme's one concentrated division (Sec 5 of the doc).
94 *
95 * All arithmetic is f64 (JS numbers) regardless of the input arrays' storage
96 * type; results are rounded to f32 only on upload. That is the doc's "CPU
97 * precompute in float64" mitigation, inherited for free.
98 */
99export interface FluxMetricFields {
100 p1: Float64Array;
101 p2: Float64Array;
102 q2: Float64Array;
103 r: Float64Array;
104}
106/**
107 * sXt* are the Cartesian components of sin(theta)*X_theta, Xp* those of
108 * X_phi, all grid space, npts each. No sin(theta) input is needed: every
109 * division the scheme performs is by D, which the sin-weighted inputs build
110 * directly.
111 */
112export function computeFluxMetric(
113 npts: number,
114 sXtx: ArrayLike<number>,
115 sXty: ArrayLike<number>,
116 sXtz: ArrayLike<number>,
117 Xpx: ArrayLike<number>,
118 Xpy: ArrayLike<number>,
119 Xpz: ArrayLike<number>,
120): FluxMetricFields {
121 const p1 = new Float64Array(npts);
122 const p2 = new Float64Array(npts);
123 const q2 = new Float64Array(npts);
124 const r = new Float64Array(npts);
126 for (let i = 0; i < npts; i++) {
127 const xt = sXtx[i];
128 const yt = sXty[i];
129 const zt = sXtz[i];
130 const xp = Xpx[i];
131 const yp = Xpy[i];
132 const zp = Xpz[i];
134 const gtt = xt * xt + yt * yt + zt * zt; // sin^2 g_tt
135 const gtp = xt * xp + yt * yp + zt * zp; // sin g_tp
136 const gpp = xp * xp + yp * yp + zp * zp; // g_pp
137 const D = Math.sqrt(gtt * gpp - gtp * gtp); // J sin^2(theta)
139 p1[i] = gpp / D;
140 p2[i] = -gtp / D;
141 q2[i] = gtt / D;
142 r[i] = 1 / D;
143 }
145 return { p1, p2, q2, r };
146}