1function prob = build_problem(a, k, p, d, wantGrid, nearBoundary)
2% BUILD_PROBLEM Assemble the problem struct for laplace-dirichlet-2d.
3%
4% The domain is bounded by one of two analytic curves, selected by p:
5%
6% p == 0 the star family, r(t) = 1 + a cos(k t),
7% p > 0 the rounded square, r(t) = (cos^p t + sin^p t)^(-1/p), which
8% traces the superellipse |x|^p + |y|^p = 1: four nearly
9% straight sides meeting at four corners of radius about 1.4/p.
10% a and k are unused. p must be even, which keeps the curve
11% analytic, so curveDD exists everywhere as the interface
12% requires.
13%
14% In both cases x(t) = r(t) [cos t; sin t], t in [0, 2 pi).
15%
16% nearBoundary (default 0) adds the near-boundary target set to the
17% evaluation points: four points per parameter just inside the boundary
18% along the inward normal, at eight parameters around the curve.
19%
20% The exact solution is u(x) = sum_j c_j log|x - s_j|, with three point
21% sources outside the domain: s_j is the boundary point at parameter
22% phi_j = 2 pi (j-1)/3 + 0.4 pushed a distance d along the outward unit
23% normal, with strengths c = [1.0; -0.6; 0.8].
24%
25% The solver receives only the curve (with derivatives), the Dirichlet
26% data g as a function of the boundary parameter, and the points where
27% the solution is requested. The sources exist here only to manufacture
28% the data; a submitted solver must not use knowledge of them.
29%
30% Fields of prob:
31% curve @(t) -> [x y] boundary point, t column vector
32% curveD @(t) -> [x' y'] first derivative
33% curveDD @(t) -> [x'' y''] second derivative
34% g @(t) -> g Dirichlet data at boundary parameter t
35% evalXY m x 2 points where uEval is required (289; 305
36% on an instance with corners, 321 on one
37% with the near-boundary set)
38% vizXY m x 2 grid points where uGrid is requested
39% (m = 0 when no visualization is wanted)
41if nargin < 6
42 nearBoundary = 0;
43end
45corners = p > 0;
47if corners
48 prob = struct();
49 prob.curve = @(t) sq_xy(t, p);
50 prob.curveD = @(t) sq_xyd(t, p);
51 prob.curveDD = @(t) sq_xydd(t, p);
52 radius = @(t) sq_r(t, p);
53 radiusD = @(t) sq_rd(t, p);
54 rmax = sq_r(pi/4, p);
55else
56 prob = struct();
57 prob.curve = @(t) [(1 + a*cos(k*t)).*cos(t), (1 + a*cos(k*t)).*sin(t)];
58 prob.curveD = @(t) [-a*k*sin(k*t).*cos(t) - (1 + a*cos(k*t)).*sin(t), ...
59 -a*k*sin(k*t).*sin(t) + (1 + a*cos(k*t)).*cos(t)];
60 prob.curveDD = @(t) [(-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*cos(t) + 2*a*k*sin(k*t).*sin(t), ...
61 (-a*k*k*cos(k*t) - 1 - a*cos(k*t)).*sin(t) - 2*a*k*sin(k*t).*cos(t)];
62 radius = @(t) 1 + a*cos(k*t);
63 radiusD = @(t) -a*k*sin(k*t);
64 rmax = 1 + abs(a);
65end
67% The three sources: boundary points at phi_j pushed d along the outward
68% normal. The normal comes from the radius and its derivative, so this is
69% the same formula for both families.
70phi = 2*pi*[0; 1; 2]/3 + 0.4;
71c = [1.0; -0.6; 0.8];
72rphi = radius(phi);
73drphi = radiusD(phi);
74bx = rphi.*cos(phi);
75by = rphi.*sin(phi);
76dxb = drphi.*cos(phi) - rphi.*sin(phi);
77dyb = drphi.*sin(phi) + rphi.*cos(phi);
78sp = sqrt(dxb.^2 + dyb.^2);
79sx = bx + d*(dyb./sp);
80sy = by - d*(dxb./sp);
82prob.g = @(t) laplace2d_bdata(t, radius(t), sx, sy, c);
84% Evaluation points: 32 rays, radial fractions 0.1..0.9, plus the origin
85% (289 points). On an instance with corners, four further points per
86% corner follow, just inside the boundary along the corner's diagonal
87% (305 points). On an instance carrying the near-boundary set, four
88% further points per parameter follow after those, at the same four
89% distances inside the boundary but along the inward unit normal, at
90% eight parameters around the curve (321 points). The rule must match
91% evalPoints() in src/problems/laplace2d/exact.ts.
92rho = (1:9)'/10;
93th = 2*pi*(0:31)'/32 + 0.13;
94deltas = [0.005; 0.01; 0.02; 0.05];
95if corners
96 cth = pi/4 + pi*(0:3)'/2;
97else
98 cth = zeros(0, 1);
99end
100if nearBoundary
101 ntt = 2*pi*(0:7)'/8 + 0.07;
102else
103 ntt = zeros(0, 1);
104end
105npts = numel(rho)*numel(th) + 1 + (numel(cth) + numel(ntt))*numel(deltas);
106pts = zeros(npts, 2);
107idx = 1;
108for i = 1:numel(rho)
109 for j = 1:numel(th)
110 rr = rho(i)*radius(th(j));
111 pts(idx, 1) = rr*cos(th(j));
112 pts(idx, 2) = rr*sin(th(j));
113 idx = idx + 1;
114 end
115end
116pts(idx, 1) = 0;
117pts(idx, 2) = 0;
118idx = idx + 1;
119for i = 1:numel(cth)
120 rc = radius(cth(i));
121 for j = 1:numel(deltas)
122 pts(idx, 1) = (rc - deltas(j))*cos(cth(i));
123 pts(idx, 2) = (rc - deltas(j))*sin(cth(i));
124 idx = idx + 1;
125 end
126end
127% The inward normal of the counterclockwise curve is -(y', -x')/|x'|.
128rn = radius(ntt);
129drn = radiusD(ntt);
130dxn = drn.*cos(ntt) - rn.*sin(ntt);
131dyn = drn.*sin(ntt) + rn.*cos(ntt);
132spn = sqrt(dxn.^2 + dyn.^2);
133for i = 1:numel(ntt)
134 for j = 1:numel(deltas)
135 pts(idx, 1) = rn(i)*cos(ntt(i)) - deltas(j)*dyn(i)/spn(i);
136 pts(idx, 2) = rn(i)*sin(ntt(i)) + deltas(j)*dxn(i)/spn(i);
137 idx = idx + 1;
138 end
139end
140prob.evalXY = pts;
142% Visualization grid: ngrid x ngrid points over the bounding square,
143% listed with y varying fastest (MATLAB column order). Points outside
144% the domain are included; the viewer masks them.
145if wantGrid
146 ngrid = 200;
147 R = 1.05*rmax;
148 xs = linspace(-R, R, ngrid);
149 [X, Y] = meshgrid(xs, xs);
150 prob.vizXY = [X(:), Y(:)];
151else
152 prob.vizXY = zeros(0, 2);
153end
155end
157% --- the rounded-square radius and its derivatives ---------------------
158%
159% With f(t) = cos^p t + sin^p t, the radius is r = f^(-1/p). For p in the
160% hundreds f underflows near the diagonals, so r goes through a logarithm
161% and the derivatives are written with the ratios f'/f and f''/f, which
162% stay of moderate size.
164function r = sq_r(t, p)
165c = cos(t); s = sin(t);
166r = exp(-log(c.^p + s.^p)/p);
167end
169function r1 = sq_rd(t, p)
170c = cos(t); s = sin(t);
171f = c.^p + s.^p;
172fp = p*(s.^(p-1).*c - c.^(p-1).*s);
173r1 = -(sq_r(t, p)/p).*(fp./f);
174end
176function r2 = sq_rdd(t, p)
177c = cos(t); s = sin(t);
178f = c.^p + s.^p;
179fp = p*(s.^(p-1).*c - c.^(p-1).*s);
180fpp = p*((p-1)*(s.^(p-2).*c.^2 + c.^(p-2).*s.^2) - f);
181u = fp./f;
182v = fpp./f;
183r2 = sq_r(t, p).*((1/p)*(1/p + 1)*u.^2 - v/p);
184end
186function xy = sq_xy(t, p)
187r = sq_r(t, p);
188xy = [r.*cos(t), r.*sin(t)];
189end
191function xy = sq_xyd(t, p)
192r = sq_r(t, p);
193r1 = sq_rd(t, p);
194xy = [r1.*cos(t) - r.*sin(t), r1.*sin(t) + r.*cos(t)];
195end
197function xy = sq_xydd(t, p)
198r = sq_r(t, p);
199r1 = sq_rd(t, p);
200r2 = sq_rdd(t, p);
201xy = [r2.*cos(t) - 2*r1.*sin(t) - r.*cos(t), ...
202 r2.*sin(t) + 2*r1.*cos(t) - r.*sin(t)];
203end