/ concept-collection / fastandaccurate
Sign in
concept-collection / fastandaccurate
fastandaccurate / src / solvers / chunkie-dlp / solver.m
64 lines · 2.0 KBCodeBlameHistory
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
e58e208chunkie-dlp runs in real MATLAB via the CLI; solver runtime field, matlab -batch harnessJeremy Magland 11% with corrected quadrature for targets near the boundary.
13% This solver runs in real MATLAB only: the command-line harness invokes
14% it through `matlab -batch` with chunkie on the path (fetched on first
15% use). It is not runnable in the browser.
17% n : number of chunks (16 points each).
19chnkr = chunkerfuncuni(@(t) fcurve(t, prob), n);
21% Dirichlet data at the nodes: for this problem family the curve
22% parameter is the polar angle, so it is recovered from the node
23% coordinates.
24xy = chnkr.r(:, :);
25tt = mod(atan2(xy(2, :), xy(1, :)), 2*pi);
26rhs = prob.g(tt(:));
28fkern = kernel('lap', 'd');
29sysmat = chunkermat(chnkr, fkern);
30sysmat = sysmat - 0.5*eye(chnkr.npt);
31sigma = sysmat \ rhs;
33out = struct();
34out.uEval = eval_targets(chnkr, fkern, sigma, prob.evalXY);
35if size(prob.vizXY, 1) > 0
36 out.uGrid = eval_targets(chnkr, fkern, sigma, prob.vizXY);
37else
38 out.uGrid = zeros(0, 1);
39end
41end
43function u = eval_targets(chnkr, fkern, sigma, XY)
44% Direct (unaccelerated) evaluation, in blocks to bound memory. accel is
e58e208chunkie-dlp runs in real MATLAB via the CLI; solver runtime field, matlab -batch harnessJeremy Magland 45% disabled so that the FLAM/fmm2d submodules are not required; at these
46% sizes direct evaluation is cheap anyway.
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
moveopenescclose