1function prob = build_problem(a, k, d, family, 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 depends on the data family:
7% family 0 (log-sources): u(x) = sum_j c_j log|x - s_j|, with three
8% point sources s_j outside the domain: s_j is the boundary point at
9% parameter phi_j = 2 pi (j-1)/3 + 0.4 pushed a distance d along the
10% outward unit normal, with strengths c = [1.0; -0.6; 0.8].
11% family 1 (branch-point): u(z) = Re sqrt(-(z - z0) e^{-i th0}), with
12% z0 the boundary point at parameter 0.4 pushed a distance d along
13% the outward unit normal and th0 its polar angle (see
14% laplace2d_bdata_branch.m).
15%
16% The solver receives only the curve (with derivatives), the Dirichlet
17% data g as a function of the boundary parameter, and the points where
18% the solution is requested. The singularities exist here only to
19% manufacture the data; a submitted solver must not use knowledge of
20% them.
21%
22% Fields of prob:
23% curve @(t) -> [x y] boundary point, t column vector
24% curveD @(t) -> [x' y'] first derivative
25% curveDD @(t) -> [x'' y''] second derivative
26% g @(t) -> g Dirichlet data at boundary parameter t
27% evalXY 289 x 2 points where uEval is required
28% vizXY m x 2 grid points where uGrid is requested
29% (m = 0 when no visualization is wanted)
31if family == 1
32 phi = 0.4;
33else
34 phi = 2*pi*[0; 1; 2]/3 + 0.4;
35end
36c = [1.0; -0.6; 0.8];
37rphi = 1 + a*cos(k*phi);
38bx = rphi.*cos(phi);
39by = rphi.*sin(phi);
40dxb = -a*k*sin(k*phi).*cos(phi) - rphi.*sin(phi);
41dyb = -a*k*sin(k*phi).*sin(phi) + rphi.*cos(phi);
42sp = sqrt(dxb.^2 + dyb.^2);
43sx = bx + d*(dyb./sp);
44sy = by - d*(dxb./sp);
46prob = struct();
47prob.curve = @(t) [(1 + a*cos(k*t)).*cos(t), (1 + a*cos(k*t)).*sin(t)];
48prob.curveD = @(t) [-a*k*sin(k*t).*cos(t) - (1 + a*cos(k*t)).*sin(t), ...
49 -a*k*sin(k*t).*sin(t) + (1 + a*cos(k*t)).*cos(t)];
50prob.curveDD = @(t) [(-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*cos(t) + 2*a*k*sin(k*t).*sin(t), ...
51 (-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*sin(t) - 2*a*k*sin(k*t).*cos(t)];
52if family == 1
53 th0 = atan2(sy(1), sx(1));
54 prob.g = @(t) laplace2d_bdata_branch(t, a, k, sx(1), sy(1), th0);
55else
56 prob.g = @(t) laplace2d_bdata(t, a, k, sx, sy, c);
57end
59% Evaluation points: 32 rays, radial fractions 0.1..0.9, plus the origin
60% (289 points). The rule must match evalPoints() in
61% src/problems/laplace2d/exact.ts.
62rho = (1:9)'/10;
63th = 2*pi*(0:31)'/32 + 0.13;
64pts = zeros(numel(rho)*numel(th) + 1, 2);
65idx = 1;
66for i = 1:numel(rho)
67 for j = 1:numel(th)
68 rr = rho(i)*(1 + a*cos(k*th(j)));
69 pts(idx, 1) = rr*cos(th(j));
70 pts(idx, 2) = rr*sin(th(j));
71 idx = idx + 1;
72 end
73end
74prob.evalXY = pts;
76% Visualization grid: ngrid x ngrid points over the bounding square,
77% listed with y varying fastest (MATLAB column order). Points outside
78% the domain are included; the viewer masks them.
79if wantGrid
80 ngrid = 200;
81 R = 1.05*(1 + abs(a));
82 xs = linspace(-R, R, ngrid);
83 [X, Y] = meshgrid(xs, xs);
84 prob.vizXY = [X(:), Y(:)];
85else
86 prob.vizXY = zeros(0, 2);
87end
89end