8144287WIP: added code for Laplace-Beltrami operator evaluation on smooth genus-0 surfaceOwen Melia 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}