/*
 * contour.c --- band2 の出力データで等高線を描く
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define G_DOUBLE
#include <glsc.h>
#include <matrix.h>

typedef struct {
    int node[3];
} threenodes;

void open_window(double xmin, double xmax, double ymin, double ymax);
void draw_contour(int, threenodes *, vector, vector, vector, double, double);
void draw_elements(int, threenodes *, vector, vector, int, int);
static void celldraw(double, double, double,
		     double, double, double,
		     double, double, double,
		     double);

double  f(double, double);
void    errorexit(char *);
void input0(FILE *fp, int *nnode, int *nelmt, int *nbc, int *nband);
void input(FILE *fp, int nnode, int nelmt, int nbc, int nband,
           threenodes *ielmt, vector x, vector y, vector u, ivector ibc);
void find_range(int nnode, vector x, double *xmin, double *xmax);

/*      main program
 *      the finite element method for the Poisson equation
 */

int verbose = 0;

int main(int argc, char **argv)
{
    /*
     * nnode        節点の総数
     * nelmt        有限要素の総数
     * nbc          基本境界条件を課す節点の総数
     * nband        半バンド幅+1 (対角行列 1, 三重対角 2, ...)
     * x[], y[]     節点の座標
     * ielmt[][3]   各有限要素を構成する節点の番号
     * ibc[]        基本境界条件を課す節点の番号
     * am[][]       全体行列の右上半分のバンド内
     * fm[]         全体ベクトル
     *
     */

    int nnode, nelmt, nbc, nband,i,j;
    double xmin, xmax, ymin, ymax, umin, umax, dx, dy;
    vector x, y, u;
    matrix am;
    ivector ibc;
    threenodes *ielmt;
    FILE *fp;
    /* 要素の縁を描くか */
    int draw_element_edge = 0;
    /* 領域に色を塗るか */
    int fill_region = 1;

    if (argc == 2)
        fp = fopen(argv[1], "r");
    else
        fp = stdin;

    /* nnode, nelmt, nbc を入力 */
    input0(fp, &nnode, &nelmt, &nbc, &nband);

    /* 大きな記憶領域を動的に確保 */
    if ((ielmt = malloc(sizeof(threenodes) * nelmt)) == NULL)
        errorexit("要素内節点の番号表 threenodes の確保が出来ませんでした。");
    if ((ibc = new_ivector(nnode)) == NULL)
        errorexit("基本境界条件節点の番号表 ibc の確保が出来ませんでした。");
    if ((x = new_vector(nnode)) == NULL)
        errorexit("節点の x 座標の表 x の確保が出来ませんでした。");
    if ((y = new_vector(nnode)) == NULL)
        errorexit("節点の y 座標の表 y の確保が出来ませんでした。");
    if ((u = new_vector(nnode)) == NULL)
        errorexit("解 u の確保が出来ませんでした。");

    /* 入力 */
    input(fp, nnode, nelmt, nbc, nband, ielmt, x, y, u, ibc);

    find_range(nnode, x, &xmin, &xmax);
    find_range(nnode, y, &ymin, &ymax);
    find_range(nnode, u, &umin, &umax);

    open_window(xmin, xmax, ymin, ymax);

    draw_elements(nelmt, ielmt, x, y, draw_element_edge, fill_region);
    draw_contour(nelmt, ielmt, x, y, u, umin, umax);

    g_sleep(-1.0);
    g_term();

    return 0;
}

void find_range(int nnode, vector x, double *xmin, double *xmax)
{
    int i;
    *xmin = *xmax = x[0];
    for (i = 0; i < nnode; i++)
        if (x[i] > *xmax)
            *xmax = x[i];
        else if (x[i] < *xmin)
            *xmin = x[i];
}

void input0(FILE *fp, int *nnode, int *nelmt, int *nbc, int *nband)
{
    fscanf(fp, "%d %d %d %d", nnode, nelmt, nbc, nband);
}

void input(FILE *fp, int nnode, int nelmt, int nbc, int nband,
           threenodes *ielmt, vector x, vector y, vector u, ivector ibc)
{
    int i, j;
    /* input */
    for (i = 0; i < nnode; i++)
        fscanf(fp, "%lf %lf", &(x[i]), &(y[i]));
    for (i = 0; i < nelmt; i++)
        for (j = 0; j < 3; j++)
            fscanf(fp, "%d", &(ielmt[i].node[j]));

    for (i = 0; i < nbc; i++)
        fscanf(fp, "%d", &ibc[i]);

    for (i = 0; i < nnode; i++)
        fscanf(fp, "%lf", &u[i]);

    /* output of input data; */
    if (verbose) {
      printf("basic data\n");
      printf("    nnode=%4d     nelmt=%4d     nbc=%4d     nband=%4d\n",
	     nnode, nelmt, nbc, nband);
      printf("x,y-coordinates of nodes (節点の x,y 座標)\n");
      printf("    i    x(i)    y(i)   i    x(i)    y(i)\n");
      for (i = 0; i < nnode; i++) {
        printf("%4d %8.4f %8.4f", i, x[i], y[i]);
        if (i % 2 == 1) printf("\n");
      }
      printf("\nnodes of elements (各要素を構成する節点)\n",
	     "    i   i1  i2  i3\n    i    i1    i2    i3\n");
      for (i = 0; i < nelmt; i++) {
        printf("%4d %4d %4d %4d ",
	       i, ielmt[i].node[0], ielmt[i].node[1], ielmt[i].node[2]);
        if (i % 2 == 1)
	  printf("\n");
      }
      if (nbc > 0) {
        printf("\nnodes with zero Dirichlet data (基本境界条件を課す節点)\n");
        for (i = 0; i < nbc; i++)
	  printf("%5d", ibc[i]);
        if (i % 8 == 7) printf("\n");
      }
      printf("\n");
      for (i = 0; i < nnode; i++) {
        printf("%8.4f ", u[i]);
        if (i % 2 == 1) printf("\n");
      }
    }
}

void errorexit(char *msg)
{
    fprintf(stderr, msg);
    exit(1);
}


void open_window(double xmin, double xmax, double ymin, double ymax)
{
  double widx, widy, x_margin, y_margin, Lx, Ly;

  x_margin = 0.1 * (xmax - xmin);
  y_margin = 0.1 * (ymax - ymin);
  widx = xmax - xmin;
  widy = ymax - ymin;

  if (widx >= widy) {
    Lx = 150.0;
    Ly = Lx * widy / widx;
  }
  else {
    Ly = 150.0;
    Lx = Ly * widx / widy;
  }
  g_init("Contour", Lx + 10.0, Ly + 10.0);
  g_device(G_BOTH);
  g_def_scale(0,
	      xmin - x_margin, xmax + x_margin,
	      ymin - y_margin, ymax + y_margin,
	      5.0, 5.0, Lx, Ly);
  g_cls();
  g_sel_scale(0);

#ifdef NONE
  g_line_color(G_BLUE);
  g_area_color(G_YELLOW);
  g_box(xmin - x_margin, xmax + x_margin, ymin - y_margin, ymax + y_margin,
	G_YES, G_YES);
#endif
}

/* 三角形要素を描く */
void draw_elements(int num_element, threenodes *ielmt,
		   vector x, vector y,
		   int draw_element_edge, int fill_region)
{
  int i, j;
  double zx[3], zy[3];

  /* 領域の縁を描く線 */
  g_def_line(0, G_BLACK, 0, G_LINE_DOTS);
  g_sel_line(0);
  /* 塗り潰す色 */
  g_def_area(0, G_YELLOW);
  g_sel_area(0);

  for (i = 0; i < num_element; i++) {
    for (j = 0; j < 3; j++) {
      zx[j] = x[ielmt[i].node[j]];
      zy[j] = y[ielmt[i].node[j]];
    }
    if (fill_region)
      g_polygon(zx, zy, 3, G_NO, G_YES);
    if (draw_element_edge) {
      g_move(zx[2], zy[2]);
      for (j = 0; j < 3; j++)
	g_plot(zx[j], zy[j]);
    }
  }
}

int level(double h[11], double u)
{
  int low, high, m;
  low = 0; high = 10;
  // h[low] <= u <= h[high] を保つ
  while (high - low >= 2) {
    m = (low + high) / 2;
    if (u >= h[m])
      low = m;
    else if (u <= h[m])
      high = m;
  }
  return low;
}

/* 等高線を描く */
void draw_contour(int num_element, threenodes *ielmt,
		  vector X, vector Y, vector U, double umin, double umax)
{
  int i, j, p, jmin, jmax;
  double x[3], y[3], u[3], max, min;
  double h[11];

  g_def_line(2, G_RED, 0, G_LINE_SOLID);
  g_sel_line(2);

  /* 等高線のレベル (最小値から最大値まで10段に分けて) */
  for (i = 0; i <= 10; i++)
    h[i] = umin + i * (umax - umin) / 10;

  for (i = 0; i < num_element; i++) {
    /* 第i要素の節点の座標, u の値を集める */
    for (j = 0; j < 3; j++) {
      p = ielmt[i].node[j];
      x[j] = X[p]; y[j] = Y[p]; u[j] = U[p];
    }
    /* u の値の範囲を求める */
    find_range(3, u, &min, &max);
    /* u の値のレベル */
    jmin = level(h, min); jmax = level(h, max);
    /* 各レベルでの等高線 */
    for (j = jmin; j <= jmax; j++)
      celldraw(x[0], x[1], x[2], y[0], y[1], y[2], u[0], u[1], u[2], h[j]);
  }
}

/* 一つの三角形内で高さ h の等高線を描く */
static void celldraw(double x1, double x2, double x3,
		     double y1, double y2, double y3,
		     double u1, double u2, double u3,
		     double h)
{
  double px[4], py[4];
  int n;
  double r;

  n = 0;

  if ((u1-h)*(u2-h) <= 0) {
    n++;
    if (u1==h) {
      px[n] = x1;
      py[n] = y1;
    }
    else if (u2 == h) {
      px[n] = x2;
      py[n] = y2;
    }
    else {
      r = (u1 - h) / (u1 - u2);
      px[n] = (1 - r) * x1 + r * x2;
      py[n] = (1 - r) * y1 + r * y2;
    }
  }
  if ((u2-h)*(u3-h) <= 0) {
    n++;
    if (u2==h) {
      px[n] = x2;
      py[n] = y2;
    }
    else if (u3 == h) {
      px[n] = x3;
      py[n] = y3;
    }
    else {
      r = (u2 - h) / (u2 - u3);
      px[n] = (1 - r) * x2 + r * x3;
      py[n] = (1 - r) * y2 + r * y3;
    }
  }
  if ((u3-h)*(u1-h) <= 0) {
    n++;
    if (u3 == h) {
      px[n] = x3;
      py[n] = y3;
    }
    else if (u1 == h) {
      px[n] = x1;
      py[n] = y1;
    }
    else {
      r = (u3 - h) / (u3 - u1);
      px[n] = (1 - r) * x3 + r * x1;
      py[n] = (1 - r) * y3 + r * y1;
    }
  }

  if (n == 2) {
    g_move(px[1], py[1]);
    g_plot(px[2], py[2]);
  }
  else if (n == 3) {
    g_move(px[1], py[1]);
    g_plot(px[2], py[2]);
    g_plot(px[3], py[3]);
    g_plot(px[1], py[1]);
  }
}