1// Exact solution, geometry, and sampling rules for laplace-dirichlet-2d,
2// implemented independently of the MATLAB side (build_problem.m). The two
3// implementations of these simple formulas check each other: a solver can
4// only reach high accuracy if both agree.
6import { hasCorners, hasNearBoundary, type Laplace2dInstance } from "./spec";
8export interface Source {
9 x: number;
10 y: number;
11 c: number;
12}
14/** Boundary radius: r(t) = 1 + a cos(k t) for the star family, and
15 * r(t) = (cos^p t + sin^p t)^(-1/p) for the rounded square, which traces
16 * the superellipse |x|^p + |y|^p = 1. Written through a logarithm because
17 * cos^p t underflows for p in the hundreds. */
18export function boundaryR(inst: Laplace2dInstance, t: number): number {
19 if (inst.shape === "rounded-square") {
20 const p = inst.p as number;
21 return Math.exp(-Math.log(Math.cos(t) ** p + Math.sin(t) ** p) / p);
22 }
23 return 1 + inst.a * Math.cos(inst.k * t);
24}
26/** dr/dt of the boundary radius. Only the ratio f'/f enters, which keeps
27 * the rounded-square case away from the underflow in f itself. */
28export function boundaryRD(inst: Laplace2dInstance, t: number): number {
29 if (inst.shape === "rounded-square") {
30 const p = inst.p as number;
31 const c = Math.cos(t);
32 const sn = Math.sin(t);
33 const f = c ** p + sn ** p;
34 const fp = p * (sn ** (p - 1) * c - c ** (p - 1) * sn);
35 return (-boundaryR(inst, t) / p) * (fp / f);
36 }
37 return -inst.a * inst.k * Math.sin(inst.k * t);
38}
40/** d2r/dt2 of the boundary radius, written like boundaryRD through the
41 * ratios f'/f and f''/f so that the rounded-square case stays away from
42 * the underflow in f itself. */
43export function boundaryRDD(inst: Laplace2dInstance, t: number): number {
44 if (inst.shape === "rounded-square") {
45 const p = inst.p as number;
46 const c = Math.cos(t);
47 const sn = Math.sin(t);
48 const f = c ** p + sn ** p;
49 const fp = p * (sn ** (p - 1) * c - c ** (p - 1) * sn);
50 const fpp =
51 p * ((p - 1) * (sn ** (p - 2) * c ** 2 + c ** (p - 2) * sn ** 2) - f);
52 const u = fp / f;
53 const v = fpp / f;
54 return boundaryR(inst, t) * ((1 / p) * (1 / p + 1) * u ** 2 - v / p);
55 }
56 return -inst.a * inst.k * inst.k * Math.cos(inst.k * t);
57}
59/** The largest radius the boundary reaches, which sets the view extent
60 * and the visualization grid. */
61export function maxRadius(inst: Laplace2dInstance): number {
62 if (inst.shape === "rounded-square") return boundaryR(inst, Math.PI / 4);
63 return 1 + Math.abs(inst.a);
64}
66/** Angles of the corners, for instances that have them: the four
67 * diagonals of the rounded square. */
68export function cornerAngles(inst: Laplace2dInstance): number[] {
69 if (!hasCorners(inst)) return [];
70 return [0, 1, 2, 3].map((j) => Math.PI / 4 + (j * Math.PI) / 2);
71}
73/** Distances inside the boundary of the near-field evaluation points. The
74 * same four distances serve the near-corner set of an instance with
75 * corners and the near-boundary set of an instance that carries one: the
76 * smallest is a third of the corner radius at p = 100, and about half a
77 * node spacing at the resolutions the sweeps reach. */
78export const NEAR_DELTAS = [0.005, 0.01, 0.02, 0.05];
80/** Boundary parameters of the near-boundary evaluation points: eight
81 * around the curve, offset so that they fall at no special phase of a lobe
82 * pattern. Empty on an instance without the set. */
83export function nearBoundaryParams(inst: Laplace2dInstance): number[] {
84 if (!hasNearBoundary(inst)) return [];
85 return Array.from({ length: 8 }, (_, j) => (2 * Math.PI * j) / 8 + 0.07);
86}
88/** Boundary point at parameter t. */
89export function boundaryPoint(inst: Laplace2dInstance, t: number) {
90 const r = boundaryR(inst, t);
91 return { x: r * Math.cos(t), y: r * Math.sin(t) };
92}
94/** The point a distance delta inside the boundary along the inward unit
95 * normal at parameter t. Since delta stays below the smallest radius of
96 * curvature on the instances that use this, the distance from the point to
97 * the curve is exactly delta. */
98export function inwardPoint(inst: Laplace2dInstance, t: number, delta: number) {
99 const r = boundaryR(inst, t);
100 const dr = boundaryRD(inst, t);
101 const dx = dr * Math.cos(t) - r * Math.sin(t);
102 const dy = dr * Math.sin(t) + r * Math.cos(t);
103 const sp = Math.hypot(dx, dy);
104 // The outward unit normal of the counterclockwise curve is (y', -x')/|x'|.
105 return {
106 x: r * Math.cos(t) - (delta * dy) / sp,
107 y: r * Math.sin(t) + (delta * dx) / sp,
108 };
109}
111/** The three exact-solution sources: boundary points at phi_j pushed a
112 * distance d along the outward normal, strengths [1.0, -0.6, 0.8]. */
113export function sources(inst: Laplace2dInstance): Source[] {
114 const { d } = inst;
115 const strengths = [1.0, -0.6, 0.8];
116 return strengths.map((c, j) => {
117 const phi = (2 * Math.PI * j) / 3 + 0.4;
118 const r = boundaryR(inst, phi);
119 const dr = boundaryRD(inst, phi);
120 const bx = r * Math.cos(phi);
121 const by = r * Math.sin(phi);
122 const dx = dr * Math.cos(phi) - r * Math.sin(phi);
123 const dy = dr * Math.sin(phi) + r * Math.cos(phi);
124 const sp = Math.hypot(dx, dy);
125 return { x: bx + (d * dy) / sp, y: by - (d * dx) / sp, c };
126 });
127}
129/** Exact solution u(x, y) = sum_j c_j log|x - s_j|. */
130export function exactU(inst: Laplace2dInstance, x: number, y: number): number {
131 let u = 0;
132 for (const s of sources(inst)) {
133 u += s.c * 0.5 * Math.log((x - s.x) ** 2 + (y - s.y) ** 2);
134 }
135 return u;
136}
138/** The evaluation points: 32 rays, radial fractions 0.1..0.9, plus the
139 * origin (289 points); on an instance with corners a further four points
140 * per corner just inside the boundary along its diagonal (305 in all); and
141 * on an instance carrying the near-boundary set four points per parameter
142 * along the inward normal (321 in all). Order matches build_problem.m:
143 * radius outer, angle inner, origin last, then the near-corner points with
144 * the corner index outer and the distance inner, then the near-boundary
145 * points with the parameter outer and the distance inner. */
146export function evalPoints(inst: Laplace2dInstance): { x: number; y: number }[] {
147 const rho = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9];
148 const pts: { x: number; y: number }[] = [];
149 for (const r of rho) {
150 for (let j = 0; j < 32; j++) {
151 const th = (2 * Math.PI * j) / 32 + 0.13;
152 const rr = r * boundaryR(inst, th);
153 pts.push({ x: rr * Math.cos(th), y: rr * Math.sin(th) });
154 }
155 }
156 pts.push({ x: 0, y: 0 });
157 for (const th of cornerAngles(inst)) {
158 const rc = boundaryR(inst, th);
159 for (const delta of NEAR_DELTAS) {
160 pts.push({ x: (rc - delta) * Math.cos(th), y: (rc - delta) * Math.sin(th) });
161 }
162 }
163 for (const t of nearBoundaryParams(inst)) {
164 for (const delta of NEAR_DELTAS) {
165 pts.push(inwardPoint(inst, t, delta));
166 }
167 }
168 return pts;
169}
171/** Exact solution at the evaluation points. */
172export function exactAtEvalPoints(inst: Laplace2dInstance): Float64Array {
173 const pts = evalPoints(inst);
174 const u = new Float64Array(pts.length);
175 pts.forEach((p, i) => (u[i] = exactU(inst, p.x, p.y)));
176 return u;
177}
179/** The visualization grid: ngrid x ngrid points over [-R, R]^2 with
180 * R = 1.05 max_t r(t). Flat index p = ix * ngrid + iy with x = xs[ix],
181 * y = xs[iy] (y varies fastest), matching build_problem.m's meshgrid
182 * column order. */
183export const VIZ_NGRID = 200;
185export function vizGrid(inst: Laplace2dInstance) {
186 const R = 1.05 * maxRadius(inst);
187 const xs = new Float64Array(VIZ_NGRID);
188 for (let i = 0; i < VIZ_NGRID; i++) {
189 xs[i] = -R + (2 * R * i) / (VIZ_NGRID - 1);
190 }
191 return { R, ngrid: VIZ_NGRID, xs };
192}
194/** Whether (x, y) is inside the domain (used only for display masking,
195 * never for scoring, so float tie-breaks at the boundary are harmless). */
196export function insideDomain(inst: Laplace2dInstance, x: number, y: number): boolean {
197 const rr = Math.hypot(x, y);
198 const th = Math.atan2(y, x);
199 return rr < boundaryR(inst, th);
200}
202/** Relative errors of numeric values against the exact solution at the
203 * evaluation points: max and L2, both relative to the exact values. */
204export function evalErrors(
205 inst: Laplace2dInstance,
206 uNum: ArrayLike<number>
207): { relMax: number; relL2: number } {
208 const uEx = exactAtEvalPoints(inst);
209 if (uNum.length !== uEx.length) {
210 throw new Error(`expected ${uEx.length} values, got ${uNum.length}`);
211 }
212 let maxDiff = 0;
213 let maxEx = 0;
214 let sumDiff2 = 0;
215 let sumEx2 = 0;
216 for (let i = 0; i < uEx.length; i++) {
217 const diff = Math.abs(uNum[i] - uEx[i]);
218 maxDiff = Math.max(maxDiff, diff);
219 maxEx = Math.max(maxEx, Math.abs(uEx[i]));
220 sumDiff2 += diff * diff;
221 sumEx2 += uEx[i] * uEx[i];
222 }
223 return {
224 relMax: maxDiff / maxEx,
225 relL2: Math.sqrt(sumDiff2 / sumEx2),
226 };
227}