% 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;

% imposing boundary condition (unlimited domain)
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

% 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);
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);
