NiHu  2.0
helmholtz_bem_2d_test.m
1 [pot_bb, pot_path] = read_epspath('potato.eps');
2 siz = norm(diff(pot_bb));
3 Le = siz/100;
4 radiator = meshpath(pot_path, Le);
5 radiator = translate_mesh(scale_mesh(radiator, 1/siz), [.2 .2]);
6 
7 field = create_slab([1.5 1], [150 100]);
8 
9 [r_nodes, r_elem] = extract_core_mesh(radiator);
10 [f_nodes, f_elem] = extract_core_mesh(field);
11 f_elem(f_elem(:,1) == 32404,1) = 22404;
12 
13 k = .7 * min(min(mesh_kmax(radiator)), min(mesh_kmax(field)));
14 [Ls, Ms, Lf, Mf] = helmholtz_bem_2d(r_nodes, r_elem, f_nodes, f_elem, k);
15 
16 x0 = [0. 0. 0.];
17 r_cent = centnorm(radiator);
18 ps = incident('line', x0, r_cent, [], k);
19 
20 f_cent = centnorm(field);
21 pf_ana = incident('line', x0, f_cent, [], k);
22 
23 I = speye(size(Ms));
24 qs_int = Ls \ (Ms - .5 * I) * ps;
25 qs_ext = Ls \ (Ms + .5 * I) * ps;
26 
27 pf = Lf * (qs_ext - qs_int);
28 
29 figure;
30 hold on;
31 plot(x0(:,1), x0(:,2), 'k*');
32 plot_mesh(field, real(pf));
33 shading flat;
34 plot_mesh(radiator);
35 axis equal tight;
36 c1 = colorbar;
37 ylabel(c1, 'Real part of numerical solution');