% local Lax-Friedrichs flux - weak form - fact that v_0 = 0 used for simplification

% computing numerical fluxes
rho_numerical_flux = ((u_0 * (rho(vmapM) + rho(vmapP)) + rho_0 * (u(vmapM) + u(vmapP))) .* n_x(:) + (rho_0 * (v(vmapM) + v(vmapP))) .* n_y(:) + (rho(vmapM) - rho(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
u_numerical_flux = ((u_0 * (u(vmapM) + u(vmapP)) + (1 / rho_0) * (p(vmapM) + p(vmapP))) .* n_x(:) + (u(vmapM) - u(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
v_numerical_flux = ((u_0 * (v(vmapM) + v(vmapP))) .* n_x(:) + ((1 / rho_0) * (p(vmapM) + p(vmapP))) .* n_y(:) + (v(vmapM) - v(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;
p_numerical_flux = ((poisson_const * p_0 * (u(vmapM) + u(vmapP)) + u_0 * (p(vmapM) + p(vmapP))) .* n_x(:) + (poisson_const * p_0 * (v(vmapM) + v(vmapP))) .* n_y(:) + (p(vmapM) - p(vmapP)) .* ((1 - alpha_flux) * C_flux(:))) / 2;

rho_numerical_flux(mapPML_M) = rho_numerical_flux(mapPML_M) + ((u_0 * (q_rho(vmapPML_M) + q_rho(vmapPML_P)) + rho_0 * (q_u(vmapPML_M) + q_u(vmapPML_P))) .* sigma_y_PML_times_n_x(mapPML_M) + (rho_0 * (q_v(vmapPML_M) + q_v(vmapPML_P))) .* sigma_x_PML_times_n_y(mapPML_M)) / 2;
u_numerical_flux(mapPML_M) = u_numerical_flux(mapPML_M) + ((u_0 * (q_u(vmapPML_M) + q_u(vmapPML_P)) + (1 / rho_0) * (q_p(vmapPML_M) + q_p(vmapPML_P))) .* sigma_y_PML_times_n_x(mapPML_M)) / 2;
v_numerical_flux(mapPML_M) = v_numerical_flux(mapPML_M) + ((u_0 * (q_v(vmapPML_M) + q_v(vmapPML_P))) .* sigma_y_PML_times_n_x(mapPML_M) + ((1 / rho_0) * (q_p(vmapPML_P) + q_p(vmapPML_P))) .* sigma_x_PML_times_n_y(mapPML_M)) / 2;
p_numerical_flux(mapPML_M) = p_numerical_flux(mapPML_M) + ((poisson_const * p_0 * (q_u(vmapPML_M) + q_u(vmapPML_P)) + u_0 * (q_p(vmapPML_M) + q_p(vmapPML_P))) .* sigma_y_PML_times_n_x(mapPML_M) + (poisson_const * p_0 * (q_v(vmapPML_M) + q_v(vmapPML_P))) .* sigma_x_PML_times_n_y(mapPML_M)) / 2;

% closing the PML with characteristic-based BC - pro tuto implementaci je potreba dopocitat vlastni vektory
%for i = 1:length(mapB)
%    node_position_faces = mapB(i);
%    node_position_total = vmapB(i);
%    variables_vector = [rho(node_position_total); u(node_position_total); v(node_position_total); p(node_position_total)];
%    numerical_flux_face = TDplusTinv_flux_cell{node_position_faces} * variables_vector;
%    rho_numerical_flux(node_position_faces) = numerical_flux_face(1);
%    u_numerical_flux(node_position_faces) = numerical_flux_face(2);
%    v_numerical_flux(node_position_faces) = numerical_flux_face(3);
%    p_numerical_flux(node_position_faces) = numerical_flux_face(4);
%end

% closing the PML with Dirichlet BC, viz. Hu 2001
% inflow
nonPMLpart = (rho_0 * u(vmapI)) .* n_x(mapI);
PMLpart = (u_0 * q_rho(vmapI) + rho_0 * q_u(vmapI)) .* sigma_y_PML_times_n_x(mapI) + (rho_0 * q_v(vmapI)) .* sigma_x_PML_times_n_y(mapI);
rho_numerical_flux(mapI) = nonPMLpart + PMLpart;
nonPMLpart = (u_0 * u(vmapI)) .* n_x(mapI);
PMLpart = (u_0 * q_u(vmapI) + (1 / rho_0) * q_p(vmapI)) .* sigma_y_PML_times_n_x(mapI);
u_numerical_flux(mapI) = nonPMLpart + PMLpart;
%nonPMLpart = zeros(length(mapI),1);
PMLpart = (u_0 * q_v(vmapI)) .* sigma_y_PML_times_n_x(mapI) + ((1 / rho_0) * q_p(vmapI)) .* sigma_x_PML_times_n_y(mapI);
%v_numerical_flux(mapI) = nonPMLpart + PMLpart; % no nonPML part
v_numerical_flux(mapI) = PMLpart;
nonPMLpart = (poisson_const * p_0 * u(vmapI)) .* n_x(mapI);
PMLpart = (poisson_const * p_0 * q_u(vmapI) + u_0 * q_p(vmapI)) .* sigma_y_PML_times_n_x(mapI) + (poisson_const * p_0 * q_v(vmapI)) .* sigma_x_PML_times_n_y(mapI);
p_numerical_flux(mapI) = nonPMLpart + PMLpart;
% outflow
nonPMLpart = (u_0 * rho(vmapO) + rho_0 * u(vmapO)) .* n_x(mapO) + (rho_0 * v(vmapO)) .* n_y(mapO);
PMLpart = (u_0 * q_rho(vmapO) + rho_0 * q_u(vmapO)) .* sigma_y_PML_times_n_x(mapO) + (rho_0 * q_v(vmapO)) .* sigma_x_PML_times_n_y(mapO);
rho_numerical_flux(mapO) = nonPMLpart + PMLpart;
nonPMLpart = (u_0 * u(vmapO)) .* n_x(mapO);
PMLpart = (u_0 * q_u(vmapO) + (1 / rho_0) * q_p(vmapO)) .* sigma_y_PML_times_n_x(mapO);
u_numerical_flux(mapO) = nonPMLpart + PMLpart;
nonPMLpart = (u_0 * v(vmapO)) .* n_x(mapO);
PMLpart = (u_0 * q_v(vmapO)) .* sigma_y_PML_times_n_x(mapO) + ((1 / rho_0) * q_p(vmapO)) .* sigma_x_PML_times_n_y(mapO);
v_numerical_flux(mapO) = nonPMLpart + PMLpart;
nonPMLpart = (poisson_const * p_0 * u(vmapO)) .* n_x(mapO) + (poisson_const * p_0 * v(vmapO)) .* n_y(mapO);
PMLpart = (poisson_const * p_0 * q_u(vmapO) + u_0 * q_p(vmapO)) .* sigma_y_PML_times_n_x(mapO) + (poisson_const * p_0 * q_v(vmapO)) .* sigma_x_PML_times_n_y(mapO);
p_numerical_flux(mapO) = nonPMLpart + PMLpart;

% resizing numerical fluxes to appropiate arrays
rho_numerical_flux = reshape(rho_numerical_flux,N_faces * N_fp,K);
u_numerical_flux = reshape(u_numerical_flux,N_faces * N_fp,K);
v_numerical_flux = reshape(v_numerical_flux,N_faces * N_fp,K);
p_numerical_flux = reshape(p_numerical_flux,N_faces * N_fp,K);

% building the right-hand side of the ODE system
rho_x = r_x .* (D_r_weak * rho) + s_x .* (D_s_weak * rho);
rho_y = r_y .* (D_r_weak * rho) + s_y .* (D_s_weak * rho);
u_x = r_x .* (D_r_weak * u) + s_x .* (D_s_weak * u);
u_y = r_y .* (D_r_weak * u) + s_y .* (D_s_weak * u);
v_x = r_x .* (D_r_weak * v) + s_x .* (D_s_weak * v);
v_y = r_y .* (D_r_weak * v) + s_y .* (D_s_weak * v);
p_x = r_x .* (D_r_weak * p) + s_x .* (D_s_weak * p);
p_y = r_y .* (D_r_weak * p) + s_y .* (D_s_weak * p);
rho_rhs = u_0 * rho_x + rho_0 * (u_x + v_y) - LIFT * (Fscale .* rho_numerical_flux) + rho_RHS(x,y,time_RK);
u_rhs = u_0 * u_x + (1 / rho_0) * p_x - LIFT * (Fscale .* u_numerical_flux);
v_rhs = u_0 * v_x + (1 / rho_0) * p_y - LIFT * (Fscale .* v_numerical_flux);
p_rhs = poisson_const * p_0 * (u_x + v_y) + u_0 * p_x - LIFT * (Fscale .* p_numerical_flux) + p_RHS(x,y,time_RK);

q_rho_x = r_x(:,elements_PML) .* (D_r_weak * q_rho(:,elements_PML)) + s_x(:,elements_PML) .* (D_s_weak * q_rho(:,elements_PML));
q_rho_y = r_y(:,elements_PML) .* (D_r_weak * q_rho(:,elements_PML)) + s_y(:,elements_PML) .* (D_s_weak * q_rho(:,elements_PML));
q_u_x = r_x(:,elements_PML) .* (D_r_weak * q_u(:,elements_PML)) + s_x(:,elements_PML) .* (D_s_weak * q_u(:,elements_PML));
q_u_y = r_y(:,elements_PML) .* (D_r_weak * q_u(:,elements_PML)) + s_y(:,elements_PML) .* (D_s_weak * q_u(:,elements_PML));
q_v_x = r_x(:,elements_PML) .* (D_r_weak * q_v(:,elements_PML)) + s_x(:,elements_PML) .* (D_s_weak * q_v(:,elements_PML));
q_v_y = r_y(:,elements_PML) .* (D_r_weak * q_v(:,elements_PML)) + s_y(:,elements_PML) .* (D_s_weak * q_v(:,elements_PML));
q_p_x = r_x(:,elements_PML) .* (D_r_weak * q_p(:,elements_PML)) + s_x(:,elements_PML) .* (D_s_weak * q_p(:,elements_PML));
q_p_y = r_y(:,elements_PML) .* (D_r_weak * q_p(:,elements_PML)) + s_y(:,elements_PML) .* (D_s_weak * q_p(:,elements_PML));
derivatives = sigma_y_PML(:,elements_PML) .* (u_0 * q_rho_x + rho_0 * q_u_x) + sigma_x_PML(:,elements_PML) .* (rho_0 * q_v_y);
non_derivatives = sigma_x_PML_plus_sigma_y_PML(:,elements_PML) .* rho(:,elements_PML) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* q_rho(:,elements_PML);
stabilization = coef_PML * (sigma_x_PML(:,elements_PML) .* (u_0 * rho(:,elements_PML) + rho_0 * u(:,elements_PML)) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* (u_0 * q_rho(:,elements_PML) + rho_0 * q_u(:,elements_PML)));;
rho_rhs(:,elements_PML) = rho_rhs(:,elements_PML) + derivatives - non_derivatives - stabilization; 
derivatives = sigma_y_PML(:,elements_PML) .* (u_0 * q_u_x + (1 / rho_0) * q_p_x);
non_derivatives = sigma_x_PML_plus_sigma_y_PML(:,elements_PML) .* u(:,elements_PML) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* q_u(:,elements_PML);
stabilization = coef_PML * (sigma_x_PML(:,elements_PML) .* (u_0 * u(:,elements_PML) + (1 / rho_0) * p(:,elements_PML)) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* (u_0 * q_u(:,elements_PML) + (1 / rho_0) * q_p(:,elements_PML)));
u_rhs(:,elements_PML) = u_rhs(:,elements_PML) + derivatives - non_derivatives - stabilization;
derivatives = sigma_y_PML(:,elements_PML) .* (u_0 * q_v_x) + sigma_x_PML(:,elements_PML) .* ((1 / rho_0) * q_p_y);
non_derivatives = sigma_x_PML_plus_sigma_y_PML(:,elements_PML) .* v(:,elements_PML) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* q_v(:,elements_PML);
stabilization = coef_PML * (sigma_x_PML(:,elements_PML) .* (u_0 * v(:,elements_PML)) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* (u_0 * q_v(:,elements_PML)));
v_rhs(:,elements_PML) = v_rhs(:,elements_PML) + derivatives - non_derivatives - stabilization;
derivatives = sigma_y_PML(:,elements_PML) .* (poisson_const * p_0 * q_u_x + u_0 * q_p_x) + sigma_x_PML(:,elements_PML) .* (poisson_const * p_0 * q_v_y);
non_derivatives = sigma_x_PML_plus_sigma_y_PML(:,elements_PML) .* p(:,elements_PML) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* q_p(:,elements_PML);
stabilization = coef_PML * (sigma_x_PML(:,elements_PML) .* (poisson_const * p_0 * u(:,elements_PML) + u_0 * p(:,elements_PML)) + sigma_x_PML_times_sigma_y_PML(:,elements_PML) .* (poisson_const * p_0 * q_u(:,elements_PML) + u_0 * q_p(:,elements_PML)));
p_rhs(:,elements_PML) = p_rhs(:,elements_PML) + derivatives - non_derivatives - stabilization;
