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

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

// advekce 1D
mat u;
double a_advec;

void Problem()
{
  T = 10;
  a_advec = 1;
  
  filename = directoryname + "P_for_h=" + mesh_size_parameter_string + ".txt";
  P.load(filename);
  Px = P.row(0);
  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);
  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
  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_1D_advection_data/"; // tj. ze ktere slozky se maji nacitat data (uloha Advection 1D)
  
  Problem();
  
  StartUp1D(); // StartUp1D jako z MATLABu
  
  BuildBCMaps();
  
  elapsed_time = timer.toc();
  cout << "preprocessing time: " << elapsed_time << " seconds" << endl;
  
  
  timer.tic();
  
  // Timestepping
  // urceni delky kroku
  double CLF = 0.75; dt = (CLF / abs(a_advec)) * min(abs(x.row(0) - x.row(1)));
  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 u_inflow;
  
  // pocatecni podminka
  double b = 1;
  u = exp((-log(2) / pow(b,2)) * square(x));
  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;
      
      
      // u_rhs = ODEs_RHS(u,time_RK);
      flux = (a_advec * u(vmapM)) % vectorise(n_x);
      
      numerical_flux = ((a_advec * (u(vmapM) + u(vmapP))) % vectorise(n_x) + (u(vmapM) - u(vmapP)) % abs(a_advec * vectorise(n_x))) / 2;
      
      u_inflow = exp((-log(2) / pow(b,2)) * square(x_in - a_advec * time_RK));
      numerical_flux(mapI) = ((a_advec * (u(vmapI) + u_inflow)) % n_x(mapI) + (u(vmapI) - u_inflow) % abs(a_advec * n_x(mapI))) / 2;
      
      numerical_flux(mapO) = (a_advec * u(vmapO)) % n_x(mapO);
      
      flux.reshape(N_fp * N_faces,K);
	  numerical_flux.reshape(N_fp * N_faces,K);
	 
	  u_rhs = -a_advec * r_x % (D_r * u) + 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 << "computation time: " << elapsed_time << " seconds" << endl;
  /*filename = "advekce_1D_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));
  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;
}
