1function out = solver(prob, n)
2% chunkie double-layer BIE for the interior Dirichlet Laplace problem.
3%
4% The same second-kind double-layer formulation as nystrom-dlp,
5% (D - I/2) sigma = g, but discretized and solved by
6% chunkie (https://github.com/fastalgorithms/chunkie), a production
7% MATLAB toolbox for boundary integral equations in 2D: the curve is
8% panelized into n uniform 16th-order Gauss-Legendre chunks, chunkermat
9% assembles the system with high-order singular quadrature, the dense
10% system is solved directly, and chunkerkerneval evaluates the potential
11% with corrected quadrature for targets near the boundary. The package
12% is fetched by mip on first use.
13%
14% n : number of chunks (16 points each).
16mip load --install magland/magland/chunkie;
18chnkr = chunkerfuncuni(@(t) fcurve(t, prob), n);
20% Dirichlet data at the nodes: for this problem family the curve
21% parameter is the polar angle, so it is recovered from the node
22% coordinates.
23xy = chnkr.r(:, :);
24tt = mod(atan2(xy(2, :), xy(1, :)), 2*pi);
25rhs = prob.g(tt(:));
27fkern = kernel('lap', 'd');
28sysmat = chunkermat(chnkr, fkern);
29sysmat = sysmat - 0.5*eye(chnkr.npt);
30sigma = sysmat \ rhs;
32out = struct();
33out.uEval = eval_targets(chnkr, fkern, sigma, prob.evalXY);
34if size(prob.vizXY, 1) > 0
35 out.uGrid = eval_targets(chnkr, fkern, sigma, prob.vizXY);
36else
37 out.uGrid = zeros(0, 1);
38end
40end
42function u = eval_targets(chnkr, fkern, sigma, XY)
43% Direct (unaccelerated) evaluation, in blocks to bound memory. accel is
44% disabled because chunkie's FMM acceleration binds to the fmm2d
45% library, which is not available in this embedded numbl runtime; at
46% these sizes direct evaluation is cheap anyway.
47opts = struct();
48opts.accel = false;
49m = size(XY, 1);
50u = zeros(m, 1);
51B = 2000;
52for i0 = 1:B:m
53 i1 = min(i0 + B - 1, m);
54 ub = chunkerkerneval(chnkr, fkern, sigma, XY(i0:i1, :).', opts);
55 u(i0:i1) = ub(:);
56end
57end
59function [r, d, d2] = fcurve(t, prob)
60tt = t(:);
61r = prob.curve(tt).';
62d = prob.curveD(tt).';
63d2 = prob.curveDD(tt).';
64end