// This file is part of the project: Einfuehrung in die Programmierung in C++
// Copyright holders: Felix Schindler
// License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

// WARNING: this implementation is not memory efficient (due to the visualization)

#ifndef FINITE_VOLUME_HH
#define FINITE_VOLUME_HH

#include <iostream>
#include <cmath>
#include <vector>
#include <fstream>


class FunctionInterface
{
public:
  virtual double evaluate(const double& x) const = 0;
};


void finite_volume_linear_transport(const FunctionInterface& initial_values)
{
  // analytical problem definitions
  // in addition to initial_values
  const double a = 0.0;
  const double b = 1.0;
  const double T = 2;
  const double velocity = 0.5; // <- has to be positive to match numerical flux below

  // discrete problem definitions
  const std::size_t I = 100;
  const std::size_t N = 200;

  const double h = (b - a)/I;
  const double delta_t = T/N;

  // create grid
  std::cout << "creating triangulation of [" << a << ", " << b << "] with " << I << " elements... " << std::flush;
  std::vector<double> centers(I, 0.0);
  for (std::size_t i = 0; i < I; ++i) {
    centers[i] = a + (i + 0.5)*h;
  }
  std::cout << "done" << std::endl;

  // reserve memory for the initial values
  std::vector<std::vector<double>> solution;
  solution.push_back(std::vector< double >(I, 0.0));

  // project initial values
  std::cout << "projecting initial values... " << std::flush;
  // therefore, walk the grid
  for (std::size_t i = 0; i < I; ++i) {
    // get the center of ith element
    const double& center = centers[i];
    // and evaluate the analytical initial values
    solution[0][i] = initial_values.evaluate(center);
  } // walk the grid
  std::cout << "done" << std::endl;

  // do the time integration
  std::cout << "doing time integration of [0, " << T << "] with " << N << " timesteps... " << std::flush;
  std::size_t n = 0;
  double t = 0.0;
  while (t < T) {
    // create storage for the next timestep
    solution.push_back(std::vector<double>(I, 0.0));

    // handle left boundary (we are on element i == 0)
    // get the solution of the last time step
    // since we do periodic boundary conditions, the neighbor of element 0 on the left is element I - 1
    double u_iminus1_n = solution[n][I - 1];
    double u_i_n = solution[n][0];
    // evaluate the flux
    // f(u) := velocity * u
    double f_u_iminus1_n = velocity * u_iminus1_n;
    double f_u_i_n = velocity * u_i_n;
    // compute numerical fluxes (upwind)
    double g_plus = f_u_i_n;
    double g_minus = f_u_iminus1_n;
    // compute the update
    solution[n + 1][0] = solution[n][0] - (delta_t/h)*(g_plus - g_minus);

    // for all inner entities
    // walk the grid
    for (std::size_t i = 1; i < I - 1; ++i) {
      // we are on element i
      // get the solution of the last time step
      u_iminus1_n = solution[n][i - 1];
      u_i_n = solution[n][i];
      // evaluate the flux
      f_u_iminus1_n = velocity * u_iminus1_n;
      f_u_i_n = velocity * u_i_n;
      // compute numerical fluxes (upwind)
      g_plus = f_u_i_n;
      g_minus = f_u_iminus1_n;
      // compute the update
      solution[n + 1][i] = solution[n][i] - (delta_t/h)*(g_plus - g_minus);
    } // walk the grid

    // handle right boundary (we are on element i == I - 1)
    // get the solution of the last time step
    // since we do periodic boundary conditions, the neighbor of element I - 1 on the right is element 0
    u_iminus1_n = solution[n][I - 2];
    u_i_n = solution[n][I - 1];
    // evaluate the flux
    f_u_iminus1_n = velocity * u_iminus1_n;
    f_u_i_n = velocity * u_i_n;
    // compute numerical fluxes (upwind)
    g_plus = f_u_i_n;
    g_minus = f_u_iminus1_n;
    // compute the update
    solution[n + 1][I - 1] = solution[n][I - 1] - (delta_t/h)*(g_plus - g_minus);

    // proceed to the next timestep
    t += delta_t;
    ++n;
  } // do the time integration
  std::cout << "done" << std::endl;

  // visualize the results
  std::cout << "visualizing... " << std::flush;
  const std::string filename_prefix = "solution";

  // first we write the solutions
  std::ofstream data_file(filename_prefix + ".dat");
  // therefore walk the grid
  for (std::size_t i = 0; i < solution[0].size(); ++i) {
    // write the center of element i
    data_file << centers[i] << " ";
    // write all timesteps for element i
    for (std::size_t n = 0; n < solution.size(); ++n) {
      data_file << solution[n][i] << " ";
    }
    data_file << std::endl;
  } // walk the grid
  data_file.close();

  // then an animation file
  std::ofstream ani_file(filename_prefix + ".ani");
  ani_file << "plot '" << filename_prefix << ".dat'"
           << " u ($1):ii+1 w l ls 1 t sprintf(\"u_h (%gs)\",ii*" << delta_t << ")\n"
           << "ii=ii+1\n"
           << "if (ii < n) reread" << std::endl;
  ani_file.close();

  // and then a gnuplot file
  std::ofstream gnuplot_file(filename_prefix + ".gnuplot");
  gnuplot_file << "reset\n"
               << "set term gif animate\n"
               << "set output \"" << filename_prefix << ".gif\"\n"
               << "n=" << N << "\n"
               << "set xrange [" << a << ":" << b << "]\n"
               << "set yrange [0:2]\n"
               << "ii=1\n"
               << "load \"" << filename_prefix << ".ani\"" << std::endl;
  gnuplot_file.close();
  std::cout << "done" << std::endl;

  // print some information
  std::cout << "\nRun 'gnuplot " << filename_prefix << ".gnuplot' to create the animation\n"
            << "and 'xdg-open " << filename_prefix << ".gif' to visualize the solution!\n" << std::endl;

}


#endif // FINITE_VOLUME_HH