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

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

// advekce 3D
mat u;
double a_advec, b_advec, c_advec;

void Problem()
{
  T = 10;
  a_advec = 1;
  b_advec = 1;
  c_advec = 1;
  
  filename = directoryname + "P_for_h=" + mesh_size_parameter_string + ".txt";
  P.load(filename);
  Px = P.row(0);
  Py = P.row(1);
  Pz = P.row(2);
  filename = directoryname + "E_for_h=" + mesh_size_parameter_string + ".txt";
  E.load(filename);
  E = E - 1; // kvuli cislovani od 0 v C++
}

void BuildBCMaps()
{
  vec dotproducts = a_advec * n_x(mapB) + b_advec * n_y(mapB) + c_advec * n_z(mapB);
  mapI = mapB(find(dotproducts < -1e-12));
  vmapI = vmapM(mapI);
  mapO = mapB(find(dotproducts > -1e-12));
  vmapO = vmapM(mapO);
}

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
  char str_from_sprintf_buffer[64];
  string Nplus1_string;
  
  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_3D_advection_data/"; // tj. ze ktere slozky se maji nacitat data (uloha Advection 3D)
  
  Problem();
  
  StartUp3D(); // StartUp3D jako z MATLABu
  
  BuildBCMaps();
  
  elapsed_time = timer.toc();
  cout << "preprocessing time: " << elapsed_time << " seconds" << endl;
  
  
  timer.tic();
  
  // Timestepping
  // urceni delky kroku
  sprintf(str_from_sprintf_buffer, "%d", N + 1);
  Nplus1_string = string(str_from_sprintf_buffer);
  filename = "dG_general_data/GLnodes_for_n=" + Nplus1_string + ".txt";
  vec GLnodes; GLnodes.load(filename); double r_min = abs(GLnodes(1) - GLnodes(0));
  dtscale3D(); double CFL = 0.125; dt = CFL * min(dtscale) * (2.0 / 3.0) * r_min / sqrt(pow(a_advec,2) + pow(b_advec,2) + pow(c_advec,2));
  N_steps = (int)ceil(T / dt); dt = T / N_steps;
  mat u_temp = zeros<mat>(N_p,K);
  
  // pomocne promenne pro ODEs_RHS
  mat u_rhs;
  mat flux;
  mat numerical_flux;
  mat x_in = x(vmapI);
  mat y_in = y(vmapI);
  mat z_in = z(vmapI);
  mat u_inflow;
  mat u_r, u_s, u_t, u_x, u_y, u_z;
  
  // pocatecni podminka
  double b = 1;
  u = exp((-log(2) / pow(b,2)) * (square(x) + square(y) + square(z)));
  double current_time = 0;
  double time_RK;
  //for(int timestep = 1; timestep <= N_steps; timestep++)
  for(int timestep = 1; timestep <= 200; timestep++) // docasne pro snizeni vypoctove narocnosti
  {
    for(int i = 1; i <= 5; i++)
    {
      time_RK = current_time + rk4c(i - 1) * dt;
      
      
      // u_rhs = ODEs_RHS(u,time_RK);
      flux = (a_advec * u(vmapM)) % vectorise(n_x) + (b_advec * u(vmapM)) % vectorise(n_y) + (c_advec * u(vmapM)) % vectorise(n_z);
      
      numerical_flux = ((a_advec * (u(vmapM) + u(vmapP))) % vectorise(n_x) + (b_advec * (u(vmapM) + u(vmapP))) % vectorise(n_y) + (c_advec * (u(vmapM) + u(vmapP))) % vectorise(n_z) + (u(vmapM) - u(vmapP)) % abs(a_advec * vectorise(n_x) + b_advec * vectorise(n_y) + c_advec * vectorise(n_z))) / 2;
      
      u_inflow = exp((-log(2) / pow(b,2)) * (square(x_in - a_advec * time_RK) + square(y_in - b_advec * time_RK) + square(z_in - c_advec * 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;
      
      numerical_flux(mapO) = (a_advec * u(vmapO)) % n_x(mapO) + (b_advec * u(vmapO)) % n_y(mapO) + (c_advec * u(vmapO)) % n_z(mapO);
      
      flux.reshape(N_fp * N_faces,K);
	  numerical_flux.reshape(N_fp * N_faces,K);
	  
	  u_r = D_r * u; u_s = D_s * u; u_t = D_t * u; u_x = r_x % u_r + s_x % u_s + t_x % u_t; u_y = r_y % u_r + s_y % u_s + t_y % u_t; u_z = r_z % u_r + s_z % u_s + t_z % u_t; // Grad3D
	  u_rhs = -a_advec * u_x - b_advec * u_y - c_advec * u_z + LIFT * (Fscale % (flux - numerical_flux));
      
      
      u_temp = rk4a(i - 1) * u_temp + dt * u_rhs;
	  u = u + rk4b(i - 1) * u_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 = "advekce_3D_u_spoctene_pro_N=" + N_string + "_a_h=" + mesh_size_parameter_string + "_Armadillo.txt";
  u.save(filename, raw_ascii);*/
  
  //mat u_exact = exp((-log(2) / pow(b,2)) * (square(x - a_advec * T) + square(y - b_advec * T) + square(z - c_advec * T)));
  mat u_exact = exp((-log(2) / pow(b,2)) * (square(x - a_advec * (200 * dt)) + square(y - b_advec * (200 * dt)) + square(z - c_advec * (200 * dt)))); // presne reseni po 200 casovych krocich
  mat errors = abs(u - u_exact);
  double E_max = max(max(errors));
  cout << "E_max = " << E_max << endl;
  double E_avg = accu(errors) / (N_p * K);
  cout << "E_avg = " << E_avg << endl;
  
  return 0;
}
