%% Neumann interior scattering problem
% This tutorial introduces how to use the toolbox |NIHU| to solve a simple
% interior acoustic scattering problem.
%
% The problem is scattering of an incident monopole wave field from the
% rigid walls of a shoe-box shaped room.

%% Surface and field point mesh
% Generate a brick surface with given side lengths and division parameter.
L = [5 4 3];    % dimensions of the room [x y z]
Le = .2;        % lengt of boundary elements
N = ceil(L/Le); % division parameter
mesh = create_brick_boundary(L, N);

%%
% Check that the mesh is directed inward by plotting the normals with the
% toolbox function |plot_elem_normals|
figure;
plot_mesh(mesh);
plot_elem_normals(mesh);
view(3);

%%
% All the normals look outside, so the elements have to be flipped.
mesh = flip_elements(mesh);

figure;
plot_mesh(mesh);
plot_elem_normals(mesh);
view(3);

%%
% The field point mesh consists of three slabs parallel with the coordinate
% planes
C = [
    .1 .1 .5
    .9 .1 .5
    .9 .9 .5
    .1 .9 .5
    ];
field1 = create_slab(C*diag(L), N([1,2]));
field2 = create_slab(C(:,[2 3 1])*diag(L), N([3,1]));
field3 = create_slab(C(:,[3 1 2])*diag(L), N([2,3]));
field = join_meshes(field1, field2, field3);

figure;
plot_mesh(field); view(3);

%% Incident wave field
% The incident wave field is the field of a monopole source located inside
% the room.
% The incident pressure field $p_{si}$ and pressure derivative field
% $q_{si}$ is computed in all surface DOF. In the present example,
% constant boundary elements are used, and the DOF locations are the
% element centers. The
% element centers and normals are obtained from the toolbox function
% |centnorm|, and the toolbox function |incident| computes the incident
% surface pressure and velocity field. For visualization purposes, the
% incident pressure field is computed in the field points (vertices of the
% field point mesh).

%%
% The frequency is determined from the element length
ratio = 7;                  % elements per wavelength
k = .9 * 2*pi/(ratio * Le); % wave number

r0 = [.3 .6 .7] * diag(L);  % monople source location

[center, normal] = centnorm(mesh);
[psi, qsi] = incident('point', r0, center, normal, k);

points = field.Nodes(:,2:4);    % extract field points (vertices)
pfi = incident('point', r0, points, [], k);

figure;
plot_mesh(mesh, real(psi));
c = colorbar; ylabel(c, 'Real incident pressure');
view(3);

figure;
plot_mesh(field, real(pfi)); shading interp;
c = colorbar; ylabel(c, 'Real incident pressure');
view(3);


%% BEM matrices and solution
% The BEM system matrices are generated by the toolbox function |bemHG|.
% The system of equations
%
% $$  \left({\bf H} - \frac{1}{2} {\bf I} \right) {\bf p}_{ss} = {\bf G} {\bf q}_{ss} $$
%
% is solved by matrix inversion. As a rigid scatterer is considered, the
% scattered velocity field is computed from the incident velocity as
%
% $$q_{ss} = -q_{si} $$
%
[H, G] = bemHG(mesh, k, 'const');   % Generate BEM matrices
qss = -qsi;                         % Boundary cond. at a rigid boundary
pss = (H-.5*eye(size(H))) \ (G * qss);    % solve the system

%% Field point pressure and total pressure field (post processing)
% The scattered field point pressure $p_{fs}$ is computed from the
% scattered surface pressure $p_{ss}$ and velocity $q_{ss}$ as
% 
% $${\bf p}_{fs} = {\bf H} {\bf p}_{ss} - {\bf G} {\bf q}_{ss}$$
% 
% where the system matrices are computed by the toolbox function |bemHG|.
[H, G] = bemHG(mesh, k, 'const', points);
pfs = H * pss - G * qss;

%%
% The total pressure field in the field points and on the surface is
% computed as
%
% $$ p_t = p_i + p_s $$
%
pst = psi + pss;
pft = pfi + pfs;

figure;
plot_mesh(mesh, real(pst));
c = colorbar; ylabel(c, 'Real pressure');
view(3);

figure;
plot_mesh(field, real(pft)); shading interp;
c = colorbar; ylabel(c, 'Real pressure');
view(3);