save_figures = 0;

disp('Plotting computed solutions:');

N_plot = 15; % should be chosen based on how fine the grid is
N_p_plot = (N_plot * N_plot + N_plot) / 2;
r_1 = linspace(-1,1,N_plot); % nodes on the face of I on axis x
s_1 = linspace(-1,1,N_plot); % nodes on the face of I on axis y
r_plot = r_1;
s_plot = -1 * ones(1,N_plot);
for i = 2:N_plot
    r_old = r_plot;
    s_old = s_plot;
    clear r_plot s_plot;
    r_temp = r_1(1:(N_plot - (i - 1)));
    s_temp = s_1(i) * ones(1,N_plot - (i - 1));
    r_plot = [r_old, r_temp];
    s_plot = [s_old, s_temp];
end

[a_plot,b_plot] = rstoab(r_plot,s_plot);
values_mode = zeros(N_p,N_p_plot);
temp = 1;
for i = 0:N
    for j = 0:(N - i)
        values_mode(temp,:) = ModalBaseValues2D(a_plot,b_plot,i,j);
        temp = temp + 1;
    end
end  
values_node = inv(V') * values_mode;


% density
fig_rho_CS = figure('Name','Density computed solution');

for i = 1:K
    element = E(:,i); v_1_ID = element(1); v_2_ID = element(2); v_3_ID = element(3);
    x_plot = 0.5 * (-(r_plot + s_plot)' * P(1,v_1_ID) + (r_plot + 1)' * P(1,v_2_ID) + (s_plot + 1)' * P(1,v_3_ID));
    y_plot = 0.5 * (-(r_plot + s_plot)' * P(2,v_1_ID) + (r_plot + 1)' * P(2,v_2_ID) + (s_plot + 1)' * P(2,v_3_ID));
    tri = delaunay(x_plot,y_plot);
    
    coef = rho(:,i)';
    values = coef * values_node;
    
    trisurf(tri, x_plot, y_plot, values); hold on;
end

shading interp;

xlabel('x'); ylabel('y'); zlabel('\rho');

colorbar('Location','EastOutside');
view(0, 90);
grid on;
axis equal;
xlim(x_axis_limits); ylim(y_axis_limits);


% velocity u
fig_u_CS = figure('Name','Velocity u computed solution');

for i = 1:K
    element = E(:,i); v_1_ID = element(1); v_2_ID = element(2); v_3_ID = element(3);
    x_plot = 0.5 * (-(r_plot + s_plot)' * P(1,v_1_ID) + (r_plot + 1)' * P(1,v_2_ID) + (s_plot + 1)' * P(1,v_3_ID));
    y_plot = 0.5 * (-(r_plot + s_plot)' * P(2,v_1_ID) + (r_plot + 1)' * P(2,v_2_ID) + (s_plot + 1)' * P(2,v_3_ID));
    tri = delaunay(x_plot,y_plot);
    
    coef = u(:,i)';
    values = coef * values_node;
    
    trisurf(tri, x_plot, y_plot, values); hold on;
end

shading interp;

xlabel('x'); ylabel('y'); zlabel('u');

colorbar('Location','EastOutside');
view(0, 90);
grid on;
axis equal;
xlim(x_axis_limits); ylim(y_axis_limits);


% velocity v
fig_v_CS = figure('Name','Velocity v computed solution');

for i = 1:K
    element = E(:,i); v_1_ID = element(1); v_2_ID = element(2); v_3_ID = element(3);
    x_plot = 0.5 * (-(r_plot + s_plot)' * P(1,v_1_ID) + (r_plot + 1)' * P(1,v_2_ID) + (s_plot + 1)' * P(1,v_3_ID));
    y_plot = 0.5 * (-(r_plot + s_plot)' * P(2,v_1_ID) + (r_plot + 1)' * P(2,v_2_ID) + (s_plot + 1)' * P(2,v_3_ID));
    tri = delaunay(x_plot,y_plot);
    
    coef = v(:,i)';
    values = coef * values_node;
    
    trisurf(tri, x_plot, y_plot, values); hold on;
end

shading interp;

xlabel('x'); ylabel('y'); zlabel('v');

colorbar('Location','EastOutside');
view(0, 90);
grid on;
axis equal;
xlim(x_axis_limits); ylim(y_axis_limits);


% pressure
fig_p_CS = figure('Name','Pressure computed solution');

for i = 1:K
    element = E(:,i); v_1_ID = element(1); v_2_ID = element(2); v_3_ID = element(3);
    x_plot = 0.5 * (-(r_plot + s_plot)' * P(1,v_1_ID) + (r_plot + 1)' * P(1,v_2_ID) + (s_plot + 1)' * P(1,v_3_ID));
    y_plot = 0.5 * (-(r_plot + s_plot)' * P(2,v_1_ID) + (r_plot + 1)' * P(2,v_2_ID) + (s_plot + 1)' * P(2,v_3_ID));
    tri = delaunay(x_plot,y_plot);
    
    coef = p(:,i)';
    values = coef * values_node;
    
    trisurf(tri, x_plot, y_plot, values); hold on;
end

shading interp;

xlabel('x'); ylabel('y'); zlabel('p');

colorbar('Location','EastOutside');
view(0, 90);
grid on;
axis equal;
xlim(x_axis_limits); ylim(y_axis_limits);


if save_figures == 1
    mesh_size_string = sprintf('%0.2f',mesh_size_parametr);
    mesh_size_string = strrep(mesh_size_string,'.',',');
    
    filename_string = ['Bajko_Problem3_rho_solution_for_h=',mesh_size_string,',_N=',num2str(N)];
    saveas(fig_rho_CS,filename_string,'pdf');
    saveas(fig_rho_CS,filename_string,'png');
    
    filename_string = ['Bajko_Problem3_u_solution_for_h=',mesh_size_string,',_N=',num2str(N)];
    saveas(fig_u_CS,filename_string,'pdf');
    saveas(fig_u_CS,filename_string,'png');
    
    filename_string = ['Bajko_Problem3_v_solution_for_h=',mesh_size_string,',_N=',num2str(N)];
    saveas(fig_v_CS,filename_string,'pdf');
    saveas(fig_v_CS,filename_string,'png');
    
    filename_string = ['Bajko_Problem3_p_solution_for_h=',mesh_size_string,',_N=',num2str(N)];
    saveas(fig_p_CS,filename_string,'pdf');
    saveas(fig_p_CS,filename_string,'png');
end