1% Driver for the boundary integral chapter.
2%
3% The editable script supplies
4%
5% names = method_names()
6% v = solve_bie(mid, n, a, nlobe, par, px, py)
7%
8% returning the computed solution at the points (px, py) inside the curve
9% r(t) = 1 + a cos(nlobe t). n is the method's resolution: boundary nodes for
10% the Nystrom method, source points for the method of fundamental solutions.
11% par is the auxiliary parameter each method needs (the radius of the source
12% circle, for the second one).
13%
14% The Dirichlet data is the trace of lp_bexact, so the exact solution is known
15% everywhere inside and the error needs no reference computation.
16p = lp_params();
17a = p.a;
18nlobe = p.nlobe;
19names = method_names();
20nm = numel(names);
22out = struct();
23out.names = names;
25% the curve itself, for the pictures
26tt = linspace(0, 2*pi, 401)';
27C = lp_curve(a, nlobe, tt);
28out.cx = C(:,1)';
29out.cy = C(:,2)';
31if strcmp(p.mode, 'field')
32 % a Cartesian grid over the bounding box, masked to the interior
33 m = 96;
34 rmax = 1 + abs(a) + 0.05;
35 xs = linspace(-rmax, rmax, m)';
36 [X, Y] = meshgrid(xs, xs);
37 th = atan2(Y, X);
38 rr = sqrt(X.^2 + Y.^2);
39 rb = 1 + a * cos(nlobe * th);
40 inside = rr < rb - 1e-9;
41 px = X(inside);
42 py = Y(inside);
44 t0 = tic;
45 v = solve_bie(p.mid, p.n, a, nlobe, p.par, px, py);
46 out.ms = toc(t0) * 1000;
47 ve = lp_bexact(px, py);
48 out.maxerr = max(abs(v - ve));
50 U = NaN(m, m);
51 E = NaN(m, m);
52 U(inside) = v;
53 E(inside) = abs(v - ve);
54 out.x = xs';
55 out.y = xs';
56 out.U = lp_rows(U);
57 out.Err = lp_rows(E);
58 out.unknowns = p.n;
59else
60 % errors on a polar lattice: four rings, the last of them close enough to
61 % the boundary to show what close evaluation costs
62 nth = 60;
63 th = 2*pi*(0:nth-1)'/nth;
64 rb = 1 + a * cos(nlobe * th);
65 rho = [0.25; 0.5; 0.75];
66 px = [rho(1)*rb .* cos(th); rho(2)*rb .* cos(th); rho(3)*rb .* cos(th)];
67 py = [rho(1)*rb .* sin(th); rho(2)*rb .* sin(th); rho(3)*rb .* sin(th)];
68 ve = lp_bexact(px, py);
69 % a ring at 98 per cent of the radius, where a layer potential is hardest
70 % to evaluate accurately
71 qx = 0.98 * rb .* cos(th);
72 qy = 0.98 * rb .* sin(th);
73 qe = lp_bexact(qx, qy);
75 ns = p.ns(:)';
76 nn = numel(ns);
77 E = zeros(nm, nn);
78 T = zeros(nm, nn);
79 Q = zeros(nm, nn);
80 EN = zeros(nm, nn);
81 ET = drv_methodsb(1, NaN(3*nm, nn), ns, a, nlobe, p.par, px, py, ve, qx, qy, qe, nm);
82 for k = 1:nm
83 E(k, :) = ET(3*k-2, :);
84 EN(k, :) = ET(3*k-1, :);
85 T(k, :) = ET(3*k, :);
86 Q(k, :) = lp_orders(E(k, :), ns);
87 end
88 out.Enear = lp_rows(EN);
89 out.ns = ns;
90 out.E = lp_rows(E);
91 out.T = lp_rows(T);
92 out.orders = lp_rows(Q);
93end
95lp_write(out);
97function ET = drv_methodsb(k, ET, ns, a, nlobe, par, px, py, ve, qx, qy, qe, nm)
98% Recursion rather than a loop, as in the other drivers: a call made from
99% inside a loop to a routine that itself loops is compiled through numbl's slow
100% path, and these are timed solves.
101if k > nm
102 return;
103end
104vw = solve_bie(k, 16, a, nlobe, par, px(1), py(1)); % warm-up
105R = drv_resb(1, NaN(3, numel(ns)), k, ns, a, nlobe, par, px, py, ve, qx, qy, qe);
106ET(3*k-2, :) = R(1, :);
107ET(3*k-1, :) = R(2, :);
108ET(3*k, :) = R(3, :);
109ET = drv_methodsb(k+1, ET, ns, a, nlobe, par, px, py, ve, qx, qy, qe, nm);
110end
112function R = drv_resb(j, R, mid, ns, a, nlobe, par, px, py, ve, qx, qy, qe)
113if j > numel(ns)
114 return;
115end
116n = ns(j);
117t0 = tic;
118v = solve_bie(mid, n, a, nlobe, par, px, py);
119el = toc(t0);
120if el < 0.005
121 el = drv_bestb(4, el, mid, n, a, nlobe, par, px, py);
122end
123q = solve_bie(mid, n, a, nlobe, par, qx, qy);
124R(1, j) = max(abs(v - ve));
125R(2, j) = max(abs(q - qe));
126R(3, j) = el * 1000;
127R = drv_resb(j+1, R, mid, ns, a, nlobe, par, px, py, ve, qx, qy, qe);
128end
130function el = drv_bestb(k, el, mid, n, a, nlobe, par, px, py)
131if k == 0
132 return;
133end
134t0 = tic;
135v = solve_bie(mid, n, a, nlobe, par, px, py);
136e2 = toc(t0);
137if e2 < el
138 el = e2;
139end
140el = drv_bestb(k-1, el, mid, n, a, nlobe, par, px, py);
141end