5127df5fastandaccurate: PDE solver benchmarks, laplace-dirichlet-2d with MFS and Nystrom DLP solversJeremy Magland 1function prob = build_problem(a, k, d, wantGrid)
2% BUILD_PROBLEM Assemble the problem struct for laplace-dirichlet-2d.
3%
4% The domain is the star-shaped region bounded by
5% x(t) = r(t) [cos t; sin t], r(t) = 1 + a cos(k t), t in [0, 2 pi).
6% The exact solution is u(x) = sum_j c_j log|x - s_j|, with three point
7% sources s_j outside the domain: s_j is the boundary point at parameter
8% phi_j = 2 pi (j-1)/3 + 0.4 pushed a distance d along the outward unit
9% normal, with strengths c = [1.0; -0.6; 0.8].
10%
11% The solver receives only the curve (with derivatives), the Dirichlet
12% data g as a function of the boundary parameter, and the points where
13% the solution is requested. The sources exist here only to manufacture
14% the data; a submitted solver must not use knowledge of them.
15%
16% Fields of prob:
17% curve @(t) -> [x y] boundary point, t column vector
18% curveD @(t) -> [x' y'] first derivative
19% curveDD @(t) -> [x'' y''] second derivative
20% g @(t) -> g Dirichlet data at boundary parameter t
21% evalXY 65 x 2 points where uEval is required
22% vizXY m x 2 grid points where uGrid is requested
23% (m = 0 when no visualization is wanted)
25phi = 2*pi*[0; 1; 2]/3 + 0.4;
26c = [1.0; -0.6; 0.8];
27rphi = 1 + a*cos(k*phi);
28bx = rphi.*cos(phi);
29by = rphi.*sin(phi);
30dxb = -a*k*sin(k*phi).*cos(phi) - rphi.*sin(phi);
31dyb = -a*k*sin(k*phi).*sin(phi) + rphi.*cos(phi);
32sp = sqrt(dxb.^2 + dyb.^2);
33sx = bx + d*(dyb./sp);
34sy = by - d*(dxb./sp);
36prob = struct();
37prob.curve = @(t) [(1 + a*cos(k*t)).*cos(t), (1 + a*cos(k*t)).*sin(t)];
38prob.curveD = @(t) [-a*k*sin(k*t).*cos(t) - (1 + a*cos(k*t)).*sin(t), ...
39 -a*k*sin(k*t).*sin(t) + (1 + a*cos(k*t)).*cos(t)];
40prob.curveDD = @(t) [(-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*cos(t) + 2*a*k*sin(k*t).*sin(t), ...
41 (-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*sin(t) - 2*a*k*sin(k*t).*cos(t)];
42prob.g = @(t) laplace2d_bdata(t, a, k, sx, sy, c);
44% Evaluation points: 16 rays, 4 radial fractions, plus the origin.
45% The rule must match evalPoints() in src/problems/laplace2d/exact.ts.
46rho = [0.25; 0.5; 0.75; 0.9];
47th = 2*pi*(0:15)'/16 + 0.13;
48pts = zeros(numel(rho)*numel(th) + 1, 2);
49idx = 1;
50for i = 1:numel(rho)
51 for j = 1:numel(th)
52 rr = rho(i)*(1 + a*cos(k*th(j)));
53 pts(idx, 1) = rr*cos(th(j));
54 pts(idx, 2) = rr*sin(th(j));
55 idx = idx + 1;
56 end
57end
58prob.evalXY = pts;
60% Visualization grid: ngrid x ngrid points over the bounding square,
61% listed with y varying fastest (MATLAB column order). Points outside
62% the domain are included; the viewer masks them.
63if wantGrid
64 ngrid = 200;
65 R = 1.05*(1 + abs(a));
66 xs = linspace(-R, R, ngrid);
67 [X, Y] = meshgrid(xs, xs);
68 prob.vizXY = [X(:), Y(:)];
69else
70 prob.vizXY = zeros(0, 2);
71end
73end