clc; clear; close all;

Globals3D;

tic;
Problem;

N = 4; % degree of polynomial representing solution within an element (local solution representation)

StartUp3D; % computing local operators and mesh information

[mapI, vmapI, mapO, vmapO] = BuildBCMaps();

u = IC(x,y,z); % evaluating initial condition
preprocessing_time = toc;

tic;
% u = Timestepping(u,T);
u_temp = zeros(N_p,K);

% time-step choice for stable scheme
GLnodes = GaussLobattoNodes(N + 1); r_min = abs(GLnodes(2) - GLnodes(1));
dtscale = dtscale3D(); C = 0.125; dt = C * min(dtscale) * (2 / 3) * r_min / sqrt(a_advec ^ 2 + b_advec ^ 2 + c_advec ^ 2);
N_steps = ceil(T / dt); dt = T / N_steps;

x_in = x(vmapI); y_in = y(vmapI); z_in = z(vmapI); % aby se nepocitalo tale znovu
current_time = 0;
for timestep = 1:N_steps
%for timestep = 1:200
    for i = 1:5
        time_RK = current_time + rk4c(i) * dt;
        
        
        % u_rhs = ODEs_RHS(u,time_RK);
        % local Lax-Friedrichs flux - strong form
        
        % computing flux on the boundaries within the elements
        flux = (a_advec * u(vmapM)) .* n_x(:) + (b_advec * u(vmapM)) .* n_y(:) + (c_advec * u(vmapM)) .* n_z(:);
        
        % computing numerical flux
        numerical_flux = ((a_advec * (u(vmapM) + u(vmapP))) .* n_x(:) + (b_advec * (u(vmapM) + u(vmapP))) .* n_y(:) + (c_advec * (u(vmapM) + u(vmapP))) .* n_z(:) + (u(vmapM) - u(vmapP)) .* abs(a_advec * n_x(:) + b_advec * n_y(:) + c_advec * n_z(:))) / 2;
        
        % imposing inflow condition
        u_inflow = BC(x_in,y_in,z_in,time_RK);
        numerical_flux(mapI) = ((a_advec * (u(vmapI) + u_inflow)) .* n_x(mapI) + (b_advec * (u(vmapI) + u_inflow)) .* n_y(mapI) + (c_advec * (u(vmapI) + u_inflow)) .* n_z(mapI) + (u(vmapI) - u_inflow) .* abs(a_advec * n_x(mapI) + b_advec * n_y(mapI) + c_advec * n_z(mapI))) / 2;
        
        % imposing outflow condition - we extrapolate u - therefore the numerical flux takes form
        numerical_flux(mapO) = (a_advec * u(vmapO)) .* n_x(mapO) + (b_advec * u(vmapO)) .* n_y(mapO) + (c_advec * u(vmapO)) .* n_z(mapO);
        
        % resizing both normal and numerical fluxes to appropiate arrays
        flux = reshape(flux,N_fp * N_faces,K);
        numerical_flux = reshape(numerical_flux,N_fp * N_faces,K);
        
        % building the right-hand side of the ODE system
        [u_x,u_y,u_z] = Grad3D(u);
        u_rhs = -a_advec * u_x - b_advec * u_y - c_advec * u_z + LIFT * (Fscale .* (flux - numerical_flux));
        
        
        
        u_temp = rk4a(i) * u_temp + dt * u_rhs;
        u = u + rk4b(i) * u_temp;
    end
    current_time = current_time + dt;
    
    if (rem(timestep,100) == 0)
        fprintf(2,'%i of total %i timesteps processed.\n',timestep,N_steps);
    end
end
computation_time = toc;