% test_ie33.m
% Test 3D infinite elements with triangular face

clear;

R = 1;
R0 = 0.4;
nR = 15;

tol = 1e-6;

ie_P = 6;           % Polynomial order (even preferred)
ie_type = 1;        % Polynomial type (1 = Lagrange)
rho0 = 1.18;        % Equilibrium density [kg/m^3]
c = 343;            % Speed of sound [m/s]


freq = 600;         % Excitation frequency

% Create the mesh
% (1) We create the boundary of a sphere quadrant boundary of triangles
mesh = create_sphere_boundary(R0, nR);
mesh = mesh_section(mesh, [-tol -tol -inf; inf inf inf]);
mesh = drop_mesh_IDs(drop_unused_nodes(mesh));
mesh = quad2tria(mesh);
% (2) Extrude the sphere quadrant manually along the radius
R_scale = ones(size(mesh.Nodes, 1), 1) * linspace(1, R/R0, nR+1);
n_nodes = size(mesh.Nodes, 1);          % Number of nodes
n_elems = size(mesh.Elements, 1);       % Number of elements
% Create new element connections
elem = mesh.Elements(:, 5:7);
shifts = repmat((0:nR-1)*n_nodes, n_elems, 1);
elem_new = [bsxfun(@plus, repmat(elem, nR, 1), shifts(:)), ...
    bsxfun(@plus, repmat(elem, nR, 1), shifts(:)+n_nodes)];
% Replace coordinates 
mesh.Nodes = [(1 : (n_nodes*(nR+1))).', ...
    bsxfun(@times, repmat(mesh.Nodes(:, 2:4), nR+1, 1), R_scale(:))];
mesh.Elements = [(1 : n_elems*nR).', ones(n_elems*nR,1)*[36 1 1], elem_new];

% Extract the outer boundary of the mesh
bou = get_boundary(mesh);
[~, ele_idx] = mesh_select(bou, sprintf('R > %f', R-tol), 'ind', 'all');
bou.Elements = bou.Elements(ele_idx, :);
bou = drop_mesh_IDs(drop_unused_nodes(bou));
% Project infinite element mesh
ie_mesh = project_ie_mesh(bou, ie_P, 'point', [0 0 0]);
% Assign infinite element properties
ie_mesh.Properties = [1 2 ie_P ie_type 0 0];

% Join the whole model
model = join_meshes(mesh, ie_mesh);
model = drop_mesh_IDs(merge_coincident_nodes(model, tol));
% Assign material
model.Materials = [1 1 rho0 c 0 0];
%% Assemble the system matrices
fprintf('Assembling system matrices ... '); tic;
[M, K, C, ~] = model_mk(model);
fprintf('Ready. (%.2f seconds)\n', toc);

%% Define an excitation
r0 = [0 0 .2];
om = 2*pi*freq;
k = om/c;
% Calculate incident field
p_in = incident('point', r0, model.Nodes(:, 2:4), [], k);

% Select excitation and solution nodes
n_nodes = size(model.Nodes, 1);
[exc_inds, ~] = mesh_select(model, sprintf('R < %f', R0+tol), 'ind');
sol_inds = setdiff((1 : n_nodes).', exc_inds);

%% Solution
S = K + 1i*om*C - om^2*M;
S21 = S(sol_inds, exc_inds);
S22 = S(sol_inds, sol_inds);

p = zeros(n_nodes, 1);
p(exc_inds) = p_in(exc_inds);
fprintf('Solving system ... '); tic;
p(sol_inds) = -S22\(S21*p_in(exc_inds));
fprintf('Ready. (%.2f seconds)\n', toc);

% Calculate error
[cmp_inds, ~] = mesh_select(model, sprintf('R < %f', R+tol), 'ind');
err = norm(p_in(cmp_inds) - p(cmp_inds)) / norm(p_in(cmp_inds));
fprintf('Relative error: %.1f %%\n', err*100);

%% Plot the solution
% Keep non-infinite elements only
[~, keep] = mesh_select(model, sprintf('R < %f', R+tol), 'ind', 'all');
model_plot = model;
model_plot.Elements = model_plot.Elements(keep, :);

figure;
subplot(1,2,1); plot_mesh(model_plot, real(p_in));
hold on; plot3(r0(1), r0(2), r0(3), 'r*');
title('Analytical solution (Real part)');
colorbar('SouthOutside')
subplot(1,2,2); plot_mesh(model_plot, real(p))
hold on; plot3(r0(1), r0(2), r0(3), 'r*');
title('Numerical solution (Real part)');

colorbar('SouthOutside');