NiHu  2.0
laplace_bem_3d_hyper_test.m
1 %// sphere radiator with radius R = 1, divided into 6 elements along the radius.
2 radiator = create_sphere_boundary(1, 6);
3 
4 %// extract the mesh description matrices
5 [r_nodes, r_elem] = extract_core_mesh(radiator);
6 
7 %// call C++ code at wave number k \approx pi
8 [Ls, Ms, Mts, Ns] = laplace_bem_3d_hyper(r_nodes, r_elem);
9 
10 %// define the incident velocity field and analytical solutions
11 x0 = [0 0 0];
12 [r_cent, r_norm] = centnorm(radiator);
13 [ps_ana, qs_ana] = incident('point', x0, r_cent, r_norm, 0);
14 
15 %// acoustic pressure on the surface and in the field points
16 ps_conv = Ms \ (Ls * qs_ana); %// conventional
17 ps_hyper = (Ns) \ (Mts * qs_ana); %// hypersingular
18 
19 % // calculate errors
20 ps_conv_err = abs(ps_conv./ps_ana - 1); %// conventional
21 ps_hyper_err = abs(ps_hyper./ps_ana - 1); %// hypersingular
22 
23 % // plot results
24 figure;
25 
26 plot_mesh(radiator, log10(ps_conv_err));
27 text(0,1.3,0, 'Conventional', 'HorizontalAlignment', 'center');
28 plot_mesh(translate_mesh(radiator, [2.3 0 0]), log10(ps_hyper_err));
29 text(2*2.3,1.3,0, 'Hypersingular', 'HorizontalAlignment', 'center');
30 shading flat;
31 view(2);
32 c = colorbar('SouthOutside');
33 xlabel(c, 'Log 10 error of solution');
34 
35 %// display error information
36 fprintf('Conventional log10 error: % .3f\n', log10(mean(ps_conv_err)));
37 fprintf('Hypersingular log10 error: % .3f\n', log10(mean(ps_hyper_err)));