/*
 * heat1d-i-glsc-eigen-sparse.cpp --- 1次元熱方程式の同次Dirichlet境界値問題
 *   Eigen の SparseMatrix を利用して3項方程式を解く
 *  2017/3/7作成
 */

#include <iostream>
#include <iomanip>
#include <vector>
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include <Eigen/SparseCholesky>
// #define G_DOUBLE
extern "C" {
#include <glsc.h>
};

using namespace std;
using namespace Eigen;

typedef Eigen::SparseMatrix<double> SpMat;
typedef Eigen::Triplet<double> T;

int main(void)
{
  int N, m, nMax, j, i, n, nskip;
  double theta, lambda, tau, h, Tmax, dt;
  double f(double x);

  cout << "N="; cin >> N;
  m = N - 1;
  h = 1.0 / N;
  cout << "lambda="; cin >> lambda;
  tau = lambda * h * h;
  cout << "theta="; cin >> theta;
  cout << "Tmax="; cin >> Tmax;
  cout << "Δt(グラフを描く時間間隔)="; cin >> dt;
  nskip = rint(dt / tau);
  // 係数行列 a を作る
  vector<T> tripletList;
  SpMat a(m,m);
  tripletList.reserve(3 * m - 2);
  for (j = 0; j < m; j++) {
    if (j == 0) {
      tripletList.push_back(T(j,   j, 1+2*theta*lambda));
      tripletList.push_back(T(j+1, j, -theta*lambda));
    }
    else if (j == m - 1) {
      tripletList.push_back(T(j-1, j, -theta*lambda));
      tripletList.push_back(T(j,   j, 1+2*theta*lambda));
    }
    else {
      tripletList.push_back(T(j-1, j, -theta*lambda));
      tripletList.push_back(T(j,   j, 1+2*theta*lambda));
      tripletList.push_back(T(j+1, j, -theta*lambda));
    }
  }
  a.setFromTriplets(tripletList.begin(), tripletList.end());
  // 係数行列を表示する
#ifdef NONE
  for (i = 0; i < a.rows(); i++) {
    for (j = 0; j < a.cols(); j++)
      cout << setw(5) << setprecision(3) << a.coeffRef(i,j) << " ";
    cout << endl;
  }
#endif
  // Cholesky分解
  SimplicialLDLT<SparseMatrix<double> > solver;
  solver.compute(a);
  if (solver.info() != Success) {
    cerr << "分解に失敗" << endl;
    exit(1);
  }
  cout << "分解できた" << endl;

  // GLSC のウィンドウを出す
  g_init((char *)"HEAT", 170.0, 170.0);
  g_device(G_BOTH);
  g_def_scale(0,
              - 0.2, 1.2, - 0.2, 1.2,
              10.0, 10.0, 150.0, 150.0);
  g_sel_scale(0);
  // 格子点上の値を記録するベクトル
  VectorXd u(N+1), b(N-1);
  Map<VectorXd> x(u.data()+1, N-1);
  // 初期値を代入し、t=0 での解のグラフを描く
  for (i = 0; i <= N; i++)
    u[i] = f(i * h);
  g_move(0.0, u[0]);
  for (i = 1; i <= N; i++)
    g_plot(i * h, u[i]);
  // 時間に関するループ
  nMax = Tmax / tau;
  for (n = 1; n <= nMax; n++) {
    for (i = 1; i < N; i++)
      b[i-1] = (1 - 2 * (1-theta) * lambda) * u[i]
        + (1 - theta) * lambda * (u[i-1] + u[i+1]);
    x = solver.solve(b);
    if (n % nskip == 0) {
      g_move(0.0, u[0]);
      for (i = 1; i <= N; i++)
        g_plot(i * h, u[i]);
    }
  }
  // クリックされるまで待つ
  g_sleep(G_STOP);
  g_term();

  return 0;
}

double f(double x)
{
  if (x < 0.5)
    return x;
  else
    return 1 - x;
}