#include <iostream>
#include <armadillo>
#include <math.h>
#include "dG_library.h"

using namespace std;
using namespace arma;
using namespace dG_global;

// informace o siti
cube TDplusTinv_flux_cube;

// LEE 2D
mat rho, u, v, p;
double rho_0, u_0, v_0, p_0, a_0, M_0;
double poisson_const = 1.4;

void Problem()
{
  	T = 75;
  	rho_0 = 1;
  	u_0 = 0;
  	v_0 = 0;
  	p_0 = rho_0 / poisson_const;
  	a_0 = sqrt(poisson_const * p_0 / rho_0);
  	M_0 = sqrt(pow(u_0,2) + pow(v_0,2)) / a_0;
  
  	filename = directoryname + "P_for_h=" + mesh_size_parameter_string + ".txt";
  	P.load(filename);
  	Px = P.row(0);
  	Py = P.row(1);
  	filename = directoryname + "E_for_h=" + mesh_size_parameter_string + ".txt";
  	E.load(filename);
  	E = E - 1; // kvuli cislovani od 0 v C++
}

int main(int argc, char* argv[])
{
  	if (argc != 3)
  	{
    	cout << "Usage: " << argv[0] << " ORDER_OF_LOCAL_APPROXIMATION MESH_SIZE_PARAMETER" << endl;
    
    	return 1;
  	}
  
  	// pomocne a debugovaci promenne
  	wall_clock timer; // mereni delky vypoctu
  	double elapsed_time;
  
  
  	timer.tic();
  
  	// parametry diskretizace
  	N = atoi(argv[1]);
  	mesh_size_parameter = atof(argv[2]);
  	DiscretizationParamatersToStrings();
  	cout << "N = " + N_string << endl;
  	cout << "h = " + mesh_size_parameter_string << endl;
  
  	directoryname = "dG_2D_CAA_Workshop1_Category4_Problem1_data/"; // tj. ze ktere slozky se maji nacitat data
  
  	Problem();
  
  	StartUp2D(); // StartUp2D jako z MATLABu
  
  	elapsed_time = timer.toc();
  	cout << "preprocessing time: " << elapsed_time << " seconds" << endl;
  
  
  	timer.tic();
  
  	// [rho,u,v,p] = Timestepping(rho,u,v,p,T);
  	// urceni delky kroku
  	filename = "dG_general_data/GLnodes_for_n=" + to_string(N_fp) + ".txt";
  	vec GLnodes; GLnodes.load(filename); double r_min = abs(GLnodes(1) - GLnodes(0));
  	dtscale2D(); double CFL = 0.75; dt = CFL * min(dtscale) * (2.0 / 3.0) * r_min / (a_0 + M_0);
  	N_steps = (int)ceil(T / dt); dt = T / N_steps;
  
  	mat rho_temp = zeros<mat>(N_p,K);
  	mat u_temp = zeros<mat>(N_p,K);
  	mat v_temp = zeros<mat>(N_p,K);
  	mat p_temp = zeros<mat>(N_p,K);
  
  	// pomocne promenne pro ODEs_RHS
  	mat rho_rhs, u_rhs, v_rhs, p_rhs;
  	mat rho_flux, rho_numerical_flux;
  	mat u_flux, u_numerical_flux;
  	mat v_flux, v_numerical_flux;
  	mat p_flux, p_numerical_flux;
  	mat u_r, u_s, u_x, u_y;
  	mat v_r, v_s, v_x, v_y;
  	mat p_r, p_s, p_x, p_y;
  	vec variables_vector, numerical_flux_face; // pro NRBC
  	int node_position_faces, node_position_total; // pro NRBC
  	mat scalar_product, u_modified, v_modified; // pro wall BC
  
  	// pocatecni podminka
  	double b = 5;
  	rho = exp((-log(2) / pow(b,2)) * (square(x) + square(y - 25)));
  	u = zeros<mat>(N_p,K);
  	v = zeros<mat>(N_p,K);
  	p = exp((-log(2) / pow(b,2)) * (square(x) + square(y - 25)));
 	double current_time = 0;
  	double time_RK;
  	for(int timestep = 1; timestep <= N_steps; timestep++)
  	{
    	for(int i = 1; i <= 5; i++)
    	{
      		time_RK = current_time + rk4c(i - 1) * dt;
      
      
      		// [rho_rhs,u_rhs,v_rhs,p_rhs] = ODEs_RHS(rho,u,v,p);
      		// computing fluxes on the boundaries within the elements
      		rho_flux = u(vmapM) % vectorise(n_x) + v(vmapM) % vectorise(n_y);
	  		u_flux = p(vmapM) % vectorise(n_x);
	  		v_flux = p(vmapM) % vectorise(n_y);
	  		p_flux = u(vmapM) % vectorise(n_x) + v(vmapM) % vectorise(n_y);
      
      		// computing numerical fluxes
	  		rho_numerical_flux = ((u(vmapM) + u(vmapP)) % vectorise(n_x) + (v(vmapM) + v(vmapP)) % vectorise(n_y) + rho(vmapM) - rho(vmapP)) / 2;
	  		u_numerical_flux = ((p(vmapM) + p(vmapP)) % vectorise(n_x) + u(vmapM) - u(vmapP)) / 2;
	  		v_numerical_flux = ((p(vmapM) + p(vmapP)) % vectorise(n_y) + v(vmapM) - v(vmapP)) / 2;
	  		p_numerical_flux = ((u(vmapM) + u(vmapP)) % vectorise(n_x) + (v(vmapM) + v(vmapP)) % vectorise(n_y) + p(vmapM) - p(vmapP)) / 2;
      
      		// imposing wall boundary condition - used also on boundaries where NRBC should be formulated
      		scalar_product = u(vmapB) % n_x(mapB) + v(vmapB) % n_y(mapB);
	  		u_modified = u(vmapB) - scalar_product % n_x(mapB);
	  		v_modified = v(vmapB) - scalar_product % n_y(mapB);
	  		rho_numerical_flux(mapB) = u_modified % n_x(mapB) + v_modified % n_y(mapB);
	  		u_numerical_flux(mapB) = p(vmapB) % n_x(mapB);
	  		v_numerical_flux(mapB) = p(vmapB) % n_y(mapB);
	  		p_numerical_flux(mapB) = u_modified % n_x(mapB) + v_modified % n_y(mapB);
      
      		// resizing both normal and numerical fluxes to appropiate arrays
      		rho_flux.reshape(N_fp * N_faces,K);
      		u_flux.reshape(N_fp * N_faces,K);
      		v_flux.reshape(N_fp * N_faces,K);
      		p_flux.reshape(N_fp * N_faces,K);
	  		rho_numerical_flux.reshape(N_fp * N_faces,K);
	  		u_numerical_flux.reshape(N_fp * N_faces,K);
	  		v_numerical_flux.reshape(N_fp * N_faces,K);
	  		p_numerical_flux.reshape(N_fp * N_faces,K);
	  
	  		//u_r = D_r * u; u_s = D_s * u; u_x = r_x % u_r + s_x % u_s; u_y = r_y % u_r + s_y % u_s; // Grad2D - staci jen u_x
	  		u_r = D_r * u; u_s = D_s * u; u_x = r_x % u_r + s_x % u_s;
	  		//v_r = D_r * v; v_s = D_s * v; v_x = r_x % v_r + s_x % v_s; v_y = r_y % v_r + s_y % v_s; // Grad2D - staci jen v_y
	  		v_r = D_r * v; v_s = D_s * v; v_y = r_y % v_r + s_y % v_s;
	  		p_r = D_r * p; p_s = D_s * p; p_x = r_x % p_r + s_x % p_s; p_y = r_y % p_r + s_y % p_s; // Grad2D
	  		rho_rhs = -u_x - v_y + LIFT * (Fscale % (rho_flux - rho_numerical_flux));
	  		u_rhs = -p_x + LIFT * (Fscale % (u_flux - u_numerical_flux));
	  		v_rhs = -p_y + LIFT * (Fscale % (v_flux - v_numerical_flux));
	  		p_rhs = -u_x - v_y + LIFT * (Fscale % (p_flux - p_numerical_flux));
      
      
      		rho_temp = rk4a(i - 1) * rho_temp + dt * rho_rhs;
      		u_temp = rk4a(i - 1) * u_temp + dt * u_rhs;
      		v_temp = rk4a(i - 1) * v_temp + dt * v_rhs;
      		p_temp = rk4a(i - 1) * p_temp + dt * p_rhs;
      
	 		rho = rho + rk4b(i - 1) * rho_temp;
	 		u = u + rk4b(i - 1) * u_temp;
	  		v = v + rk4b(i - 1) * v_temp;
	  		p = p + rk4b(i - 1) * p_temp;
    	}
    	current_time = current_time + dt;
    
    	if(timestep % 100 == 0)
    	{
      		//cout << timestep << " of total " << N_steps << " timesteps processed." << endl;
      
      		elapsed_time = timer.toc();
      		cout << timestep << " of total " << N_steps << " timesteps processed; elapsed " << elapsed_time << " seconds so far" << endl;
    	}
  	}
  
  	elapsed_time = timer.toc();
  	cout << "computation time: " << elapsed_time << " seconds" << endl;
  	filename = "CAA_Workshop1_Category4_Problem1_p_spoctene_pro_N=" + N_string + "_a_h=" + mesh_size_parameter_string + "_Armadillo.txt";
  	p.save(filename, raw_ascii);
  
  	return 0;
}
