// Dieses Programm integriert die Wärmeleitungsgleichung mit dem aus der
// Vorlesung bekannten Crank-Nicholson-Verfahren. Für die lineare Algebra wird
// dabei die Bibliothek GLS (Gnu Scientific-Library) verwendet (linken deshalb
// mit -dislin -lgsl -lgslcblas). Die Ergebnisse werden in Form von PNG-Grafiken
// ausgegeben, aus denen sich mit Programmen wie ffmpeg Videos erzeugen lassen.
// Die Anzahl der ausgegebenen Bilder pro simulierter Zeiteinheit läßt sich über
// die Variable "fps" steuern (siehe unten). 

#include <math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>

// Die folgende Zeile muss evtl. an Eure lokale Dislin-Installation
// angepasst werden.
#include "/usr/local/dislin/dislin.h"

// Eine nützliche Definition zum Runden
#define round(x) ((x)>=0?(long)((x)+0.5):(long)((x)-0.5))

// Zeitintervall [t_0,t_0]
double t_0 = 0.0;
double t_1 = 5.0;
// Raumbereich [x_0,x_1]x[y_0,y_1]
double x_0 = 0.0;
double x_1 = 1.0;

// CFL-Zahl
double alpha = 0.2;

// Nur der Raumschritt wird explizit angegeben. Wie in der Vorlesung besprochen,
// wird der Zeitschritt weiter unten aus alpha berechnet.
double dx = 0.01;

// Die Temperaturen am linken und rechten Rand des Stabes. Auskommentiert sind
// die Randwerte des zweiten Beispiels aus der Vorlesung. 
double T_l = 0.0; //2.0; 
double T_r = 0.0; //0.5;

// Die Diffusionskonstante
double D = 0.1;

// fps ist die Anzahl der Bilder, die pro simulierter Sekunde ausgegeben werden
// sollen 
int fps = 20;

// Die Funktion startwert wird zu Beginn des Programms verwendet, um die
// Anfangsbedingung festzulegen. 
double startwert(int i) { double x = x_0 + i * dx;

  // Beispiel Nr. 1 aus der Vorlesung
  return 8.0 * x * (1.0-x);

  // Beispiel Nr. 2 aus der Vorlesung
  // return 2.0 - 1.5 * x + sin(M_PI*x);
}

int main () {
 // Der Zeitschritt wird aus dx, D und der CFL-Zahl berechnet
  double dt = alpha * pow(dx, 2) / D;
  
  // Die Anzahl der räumlichen und zeitlichen Iterationsschritte wird
  // aus der Größe der entsprechenden Intervalle und den Schrittweiten
  // berechnet.
  int x_steps = round( (x_1 - x_0) / dx); 
  int t_steps = round( (t_1 - t_0) / dt);

  // Je nach Wahl der globalen Variablen fps und dt wird die Zahl der
  // zeitlichen Iterationen berechnet, die zwischen zwei Ausgaben
  // durchgeführt werden sollen.
  int steps_pro_frame = round(1.0 / fps / dt);
  if ( !steps_pro_frame ) {
    steps_pro_frame = 1;
  }
  // Die Variablen fname und titel werden für die Ausgabe mit Dislin
  // benötigt (siehe unten).
  char fname[80];
  char titel[80];

  // Für die Zeitintegration wird ein lineares Gleichungssystem der Form A x = B
  // y + 2*alpha (T_l,0,...,0,T_r) aufgestellt. Die dafür notwendigen Vektoren
  // und Gleichungen werden im Folgenden definiert. Da der Rechte und der linke
  // Rand jeweils einen festen Wert haben, haben die Vektoren die Dimension
  // x_steps-1 und die Matrizen entsprechend die Dimension
  // (x_steps-1)x(x_steps-1).
  gsl_vector *u = gsl_vector_alloc(x_steps - 1); 
  gsl_vector *temp = gsl_vector_alloc(x_steps - 1);
  gsl_matrix *A = gsl_matrix_alloc(x_steps - 1, x_steps - 1);
  gsl_matrix *B = gsl_matrix_alloc(x_steps - 1, x_steps - 1);

  // Die folgenden beiden Größen werden von der GSL intern benötigt, um die
  // linearen Gleichungen mittels LU-Zerlegung zu lösen.
  gsl_permutation *p;
  int signum;
 
  int i = 0; // i ist der x-Index
  int k = 0; // k ist der Zeitindex

  // Das Feld wird entsprechend der Anfangsbedingung mit startwert()
  // initialisiert.
  for (i = 0; i <= x_steps-2; i++) {
      gsl_vector_set(u, i, startwert(i+1));
  }

  // Erstelle die Matrix A (siehe Vorlesung)
   for ( i=0; i<=x_steps-2; i++ ) {
    gsl_matrix_set(A, i, i, 2.0 * (1+alpha));

    if (i != 0) {
      gsl_matrix_set(A, i, i-1, -alpha);
    } 
    
    if (i != x_steps-2) {
      gsl_matrix_set(A, i, i+1, -alpha);
    } 
  }

  // Erstelle die Matrix B (siehe Vorlesung)
  gsl_matrix_alloc(x_steps - 1, x_steps - 1);
  for ( i=0; i<=x_steps-2; i++ ) {
    gsl_matrix_set(B, i, i, 2.0 * (1-alpha));

    if (i != 0) {
      gsl_matrix_set(B, i, i-1, alpha);
    } 
    
    if (i != x_steps-2) {
      gsl_matrix_set(B, i, i+1, alpha);
    } 
  }

  // Als Vorbereitung auf die wiederholte Lösung des linearen Gleichungssystems
  // wird die LU-Zerlegung von A durchgeführt.
  p = gsl_permutation_alloc(x_steps - 1);
  gsl_linalg_LU_decomp(A, p, &signum); 

  // Durchlaufe jetzt alle Zeitschritte
  for ( k = 0; k <= t_steps; k++ ) {

    //---{ Ausgabe mit Dislin }-----------------------------------------
    if ( k % steps_pro_frame == 0) {
      scrmod("REVERSE");
      metafl("PNG");
      filmod("DELETE");
      sprintf(fname, "HCN_frame%d.png", round(k / steps_pro_frame) );
      setfil(fname);
      disini();
      texmod("ON");
      
      sprintf(titel, "%d Stuetzstellen, $t=%f$",x_steps + 1, t_0 + (k + 1) * dt);
      titlin(titel, 4);
    
      graf(x_0, x_1, x_0, (x_1 - x_0)/10.0, 0.0, 3.0, 0.0, 0.5);
      title();

      hsymbl(6);
      rlsymb(21, x_0, T_l);
      rlsymb(21, x_1, T_r);
      for (i=0; i<=x_steps-2; i++) {
        rlsymb(21, x_0 + (i+1) * dx, gsl_vector_get(u, i));
      }

      disfin();
   }
    //------------------------------------------------------------------


    //---{ Der eigentliche Zeitschritt }--------------------------------
    // Zunächst wird u^n mit der Matrix B multipliziert...
    gsl_blas_dgemv(CblasNoTrans, 1.0, B, u, 0.0, temp);

    // ...dann werden die Randbedingungen berücksichtigt...
    gsl_vector_set(temp, 0, gsl_vector_get(temp,0) + 2.0 * alpha * T_l);
    gsl_vector_set(temp, x_steps-2, gsl_vector_get(temp,x_steps-2) + 2.0 * alpha * T_r);

    // ...und schließlich wird das resultierende LGS gelöst.
    gsl_linalg_LU_solve(A, p, temp, u);
    //------------------------------------------------------------------

  }
  
  // Zuletzt wird der reservierte Speicher wieder freigegeben.
  gsl_vector_free(u);
  gsl_vector_free(temp);
  gsl_permutation_free(p);
  gsl_matrix_free(A);
  gsl_matrix_free(B);

  return 1;
}