1// The problem as a TypeScript object, for solvers that do not run through
2// numbl. It carries exactly what build_problem.m hands a MATLAB solver —
3// the curve with its derivatives, the Dirichlet data as a function of the
4// boundary parameter, and the points where values are required — as plain
5// functions and arrays. The specification
6// (docs/problems/laplace-dirichlet-2d.md) states the interface once and
7// this is its second form; a solver written against it, like the WebGPU
8// MFS, sees the same information as a MATLAB one and no more. In
9// particular the sources of the exact solution are used here only to
10// manufacture g, and are not reachable from the returned object.
12import { type Laplace2dInstance } from "./spec";
13import {
14 boundaryR,
15 boundaryRD,
16 boundaryRDD,
17 evalPoints,
18 sources,
19 vizGrid,
20 VIZ_NGRID,
21} from "./exact";
23export interface Vec2 {
24 x: number;
25 y: number;
26}
28export interface Laplace2dProblem {
29 /** Boundary point at parameter t. */
30 curve(t: number): Vec2;
31 /** First derivative of the curve with respect to t. */
32 curveD(t: number): Vec2;
33 /** Second derivative. */
34 curveDD(t: number): Vec2;
35 /** Dirichlet data at boundary parameter t. */
36 g(t: number): number;
37 /** The evaluation points: nEval rows of (x, y), interleaved. */
38 evalXY: Float64Array;
39 nEval: number;
40 /** The visualization grid points, interleaved, empty when not wanted. */
41 vizXY: Float64Array;
42 nViz: number;
43}
45/** Point and derivatives of x(t) = r(t) (cos t, sin t). */
46function curveAt(inst: Laplace2dInstance, t: number, order: 0 | 1 | 2): Vec2 {
47 const c = Math.cos(t);
48 const s = Math.sin(t);
49 const r = boundaryR(inst, t);
50 if (order === 0) return { x: r * c, y: r * s };
51 const r1 = boundaryRD(inst, t);
52 if (order === 1) return { x: r1 * c - r * s, y: r1 * s + r * c };
53 const r2 = boundaryRDD(inst, t);
54 return {
55 x: r2 * c - 2 * r1 * s - r * c,
56 y: r2 * s + 2 * r1 * c - r * s,
57 };
58}
60export function buildProblem(
61 inst: Laplace2dInstance,
62 wantGrid = false
63): Laplace2dProblem {
64 // The three sources exist only to manufacture g, exactly as in
65 // build_problem.m, and stay in this closure.
66 const src = sources(inst);
67 const pts = evalPoints(inst);
68 const evalXY = new Float64Array(2 * pts.length);
69 pts.forEach((p, i) => {
70 evalXY[2 * i] = p.x;
71 evalXY[2 * i + 1] = p.y;
72 });
74 let vizXY = new Float64Array(0);
75 if (wantGrid) {
76 const { xs } = vizGrid(inst);
77 vizXY = new Float64Array(2 * VIZ_NGRID * VIZ_NGRID);
78 // Flat index p = ix * ngrid + iy, y varying fastest, matching
79 // build_problem.m's meshgrid column order.
80 let k = 0;
81 for (let ix = 0; ix < VIZ_NGRID; ix++) {
82 for (let iy = 0; iy < VIZ_NGRID; iy++) {
83 vizXY[k++] = xs[ix];
84 vizXY[k++] = xs[iy];
85 }
86 }
87 }
89 return {
90 curve: (t) => curveAt(inst, t, 0),
91 curveD: (t) => curveAt(inst, t, 1),
92 curveDD: (t) => curveAt(inst, t, 2),
93 g: (t) => {
94 const p = curveAt(inst, t, 0);
95 let u = 0;
96 for (const s of src) {
97 u += s.c * 0.5 * Math.log((p.x - s.x) ** 2 + (p.y - s.y) ** 2);
98 }
99 return u;
100 },
101 evalXY,
102 nEval: pts.length,
103 vizXY,
104 nViz: vizXY.length / 2,
105 };
106}