/* pde3_gl.c  just checking second order approximation      */
/*      z(x,y) = x^3  + y^3  + 2x^2y + 3xy^2 + 4x + 5y + 6  */
/*       dz/dx = 3x^2 + 0    + 4xy   + 3y^2  + 4  + 0  + 0  */
/*   d^2z/dx^2 = 6x   + 0    + 4y    + 0     + 0  + 0  + 0  */
/*       dz/dy = 0    + 3y^2 + 2x^2  + 6xy   + 0  + 5  + 0  */
/*   d^2z/dy^2 = 0    + 6y   + 0     + 6x    + 0  + 0  + 0  */
/*                                                          */
/* solve d^2z/dx^2 + d^2z/dy^2 = 12x + 10y = c(x,y)         */
/* by second order method on square -1,-1 to 1,1            */

/* second order numerical difference                        */
/* d^2z/dx^2 + d^2z/dy^2 =                                  */
/* 1/h^2(z(-1h,0)+z(1h,0)+z(0,-1h)+z(0,1h)-4z(0,0)) +O(h^2) */
/*  = c(x,y)                                                */
/* solve for new z(0,0) =                                   */
/* ((z(-1h,0)+z(1h,0)+z(0,-1h)+z(0,1h)) - c(x,y)*h^2)/4     */

/* Mouse buttons control axis of rotation      */
/* arrow keys control axis of rotation, also   */
/* "f" and "s" make rotation faster or slower  */
/* "W" and "w" make wireframe or solid         */

#include <stdio.h>
#include <stdlib.h>
#include <GL/glut.h>

#undef  abs
#define abs(a) ((a)<0.0?(-(a)):(a))
#undef  min
#define min(a,b) ((a)<(b)?(a):(b))

static int nx = 41;
static int ny = 41;
static double zs[100][100]; /* solution being interated       */
static double zt[100][100]; /* temporary for partical results */
                 /* nx-2 by ny-2 internal, non boundary cells */
static double xmin = -2.0;
static double ymin = -2.0;
static double xmax =  4.0;
static double ymax =  4.0;
static double zmax;
static double zmin;
static double h;             /* x  and  y step = (xmax-xmin)/(n-1)     */
static double diff = 1.0e6;  /* sum abs difference form last ineration */
static double error;         /* sum of absolute exact-computed         */
static double avgerror;      /* average error                          */
static double maxerror;      /* maximum cell error                     */
static int iter = 0;         /* iteration counter                      */
static int max_iter = 3500;  /* maximum iterations                     */

/* graphics values */
static GLfloat theta[] = {60.0,90.0,0.0};
static GLint axis = 1;
static GLint windoww = 700;
static GLint windowh = 700;
static double rotation = 0.0;  /* 'f' or 's' */
static int wireframe = 1;      /* 'W' or 'w' */
static double scalex = 1.0;    /* scale x to -1..1 */
static double scaley = 1.0;    /* scale y to -1..1 */
static double scalez = 1.0;    /* scale z to -1..1 */
static double offsetx = 0.0;   /* offset to center */
static double offsety = 0.0;   /* offset to center */
static double offsetz = 0.0;   /* offset to center */
static int run = 0;            /* run/pause        */
static int step = 0;           /* step             */
static int xm, ym;             /* global mouse x,y */

/* no function prototypes used here because the functions are defined */
/* before they are used */


static double z(double x, double y) /* for boundary and check */
{
  return x*x*x + y*y*y + 2.0*x*x*y + 3.0*x*y*y + 4.0*x + 5.0*y +6.0;
}

static double c(double x, double y) /* from solve */
{
  return 12.0*x + 10.0*y;
}

static double z00(int i, int j) /* new iterated value */
{
  return ( zs[i-1][j] + zs[i+1][j] + zs[i][j-1] + zs[i][j+1] -
           c(xmin+i*h,ymin+j*h)*h*h )/4.0;
}

static void init_boundary()
{
  int i, j;
  
  zmin = z(xmin, ymin); /* save zmin and zmax for plotting */
  zmax = zmin;
  /* set boundary on four sides */
  for(i=0; i<nx; i++)
  {
    zs[i][0]    = z(xmin+i*h, ymin);
    zs[i][ny-1] = z(xmin+i*h, ymax);
    if(zs[i][0]<zmin) zmin = zs[i][0];
    if(zs[i][0]>zmax) zmax = zs[i][0];
    if(zs[i][ny-1]<zmin) zmin = zs[i][ny-1];
    if(zs[i][ny-1]>zmax) zmax = zs[i][ny-1];
  }
  for(j=0; j<ny; j++)
  {
    zs[0][j]    = z(xmin, ymin+j*h);
    zs[nx-1][j] = z(xmax, ymin+j*h);
    if(zs[0][j]<zmin) zmin = zs[0][j];
    if(zs[0][j]>zmax) zmax = zs[0][j];
    if(zs[nx-1][j]<zmin) zmin = zs[0][j];
    if(zs[nx-1][j]>zmax) zmax = zs[nx-1][j];
  }
  /* zero internal cells */
  for(i=1; i<nx-1; i++)
    for(j=1; j<ny-1; j++)
      zs[i][j] = 0.0;
} /* end init_boundary */

static void iterate()
{
  int i, j;
  
  diff = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      zt[i][j] = z00(i,j); 
      diff = diff + abs(zt[i][j]-zs[i][j]);
    }
  }
  /* copy temp to solution */
  for(i=1; i<nx-1; i++)
    for(j=1; j<ny-1; j++)
      zs[i][j] = zt[i][j];
} /* end iterate */

static void check() /* typically not known, instrumentation */
{
  int i, j;
  double zexact, err;
  
  error = 0.0;
  maxerror = 0.0;
  for(i=1; i<nx-1; i++)
  {
    for(j=1; j<ny-1; j++)
    {
      zexact = z(xmin+i*h,ymin+j*h);
      err = abs(zs[i][j]-zexact);
      error = error + err;
      if(err>maxerror) maxerror=err;
    }
  }
} /* end check */

static void printzs()
{
  int i, j;
  
  printf("\n computed solution\n");
  for(i=0; i<min(nx,11); i++)
  {
    for(j=0; j<min(ny,11); j++)
    {
      printf("%6.3f ", zs[i][j]);
    }
    printf("\n");
  }
  printf("\n");
} /* end printzs */

static void printexact()
{
  int i, j;
  
  printf("\n exact solution\n");
  for(i=0; i<min(nx,11); i++)
  {
    for(j=0; j<min(ny,11); j++)
    {
      printf("%6.3f ", z(xmin+i*h,ymin+j*h));
    }
    printf("\n");
  }
  printf("\n");
} /* end printexact */

static void printdiff()
{
  int i, j;
  
  printf("\n difference exact-computed solution\n");
  printf("\n"); 
  for(i=1; i<min(nx-1,10); i++)
  {
    printf("       ");
    for(j=1; j<min(ny-1,10); j++)
    {
      printf("%7.4f", z(xmin+i*h,ymin+j*h)-zs[i][j]);
    }
    printf("\n");
  }
  printf("\n");
} /* end printdiff */

static void draw_pde(void)
{
  int i, j;
  double x1, x2, y1, y2;
  
  /* draw surface of solution */
  if(!wireframe)
  {
    glColor3f(1.0, 0.0, 0.0);
    for(i=0; i<nx-1; i++)
    {
      for(j=0; j<ny-1; j++)
      {
        x1 = scalex*((xmin+i*h)+offsetx);
        x2 = scalex*((xmin+(i+1)*h)+offsetx);
        y1 = scaley*((ymin+j*h)+offsety);
        y2 = scaley*((ymin+(j+1)*h)+offsety);
        glBegin(GL_POLYGON);
          glVertex3f(x1, y1, scalez*(zs[i][j]+offsetz));
          glVertex3f(x2, y1, scalez*(zs[i+1][j]+offsetz));
          glVertex3f(x2, y2, scalez*(zs[i+1][j+1]+offsetz));
          glVertex3f(x1, y2, scalez*(zs[i][j+1]+offsetz));
        glEnd();
      }
    }
  }
  else /* wireframe, lines */
  {
    glLineWidth(1.0);
    glColor3f(1.0, 0.5, 0.0); /* exact, orange */
    glBegin(GL_LINES);
      for(i=0; i<nx; i++)
      {
        for(j=0; j<ny; j++)
        {
          x1 = scalex*((xmin+i*h)+offsetx);
          y1 = scaley*((ymin+j*h)+offsety);
          if(i+1<nx)
          {
            x2 = scalex*((xmin+(i+1)*h)+offsetx);
            glVertex3f(x1, y1, scalez*(z(xmin+i*h,ymin+j*h)+offsetz));
            glVertex3f(x2, y1, scalez*(z(xmin+(i+1)*h,ymin+j*h)+offsetz));
          }
          if(j+1<ny)
          {
            y2 = scaley*((ymin+(j+1)*h)+offsety);
            glVertex3f(x1, y1, scalez*(z(xmin+i*h,ymin+j*h)+offsetz));
            glVertex3f(x1, y2, scalez*(z(xmin+i*h,ymin+(j+1)*h)+offsetz));
          }
        }
      }
    glEnd();
    glLineWidth(1.0);
    glColor3f(0.0, 0.0, 0.0); /* computer solution, black */
    glBegin(GL_LINES);
      for(i=0; i<nx; i++)
      {
        for(j=0; j<ny; j++)
        {
          x1 = scalex*((xmin+i*h)+offsetx);
          y1 = scaley*((ymin+j*h)+offsety);
          if(i+1<nx)
          {
            x2 = scalex*((xmin+(i+1)*h)+offsetx);
            glVertex3f(x1, y1, scalez*(zs[i][j]+offsetz));
            glVertex3f(x2, y1, scalez*(zs[i+1][j]+offsetz));
          }
          if(j+1<ny)
          {
            y2 = scaley*((ymin+(j+1)*h)+offsety);
            glVertex3f(x1, y1, scalez*(zs[i][j]+offsetz));
            glVertex3f(x1, y2, scalez*(zs[i][j+1]+offsetz));
          }
        }
      }
    glEnd();
  }
} /* end draw_pde */

static void printstring(void *font, char *string)
{
  int len, i;

  len = (int) strlen(string);
  for (i = 0; i < len; i++) glutBitmapCharacter(font, string[i]);
} /* end printstring */

static void button(int x1, int y1, int x2, int y2, char * text)
{
  glColor3f(0.0, 0.0, 0.0);
  glRasterPos2i(x1+7, y1+7);
  printstring(GLUT_BITMAP_HELVETICA_18, text);
  glColor3f(1.0, 0.5, 0.0);
  glBegin(GL_POLYGON);
    glVertex2i(x1, y1);
    glVertex2i(x2, y1);
    glVertex2i(x2, y2);
    glVertex2i(x1, y2);
  glEnd();
}

static void printhelp(void)
{
  char iter_val[] = "iteration  =      ";
  char err_val[]  = "max error  =          ";
  char xrot_val[] = "X rotation =       degrees";
  char yrot_val[] = "Y rotation =       degrees";
  char zrot_val[] = "Z rotation =       degrees";
  char xmin_val[] = "xmin  =          ";
  char xmax_val[] = "xmax  =          ";
  char ymin_val[] = "ymin  =          ";
  char ymax_val[] = "ymax  =          ";
  char zmin_val[] = "zmin  =          ";
  char zmax_val[] = "zmax  =          ";

  sprintf(&iter_val[13], "%4d", iter);
  sprintf(&err_val[13],  "%8.3f", maxerror);
  sprintf(&xmin_val[8],  "%8.3f", xmin);
  sprintf(&xmax_val[8],  "%8.3f", xmax);
  sprintf(&ymin_val[8],  "%8.3f", ymin);
  sprintf(&ymax_val[8],  "%8.3f", ymax);
  sprintf(&zmin_val[8],  "%8.3f", zmin);
  sprintf(&zmax_val[8],  "%8.3f", zmax);
  sprintf(&xrot_val[13], "%4.0f", theta[0]);
  xrot_val[17] = ' ';
  sprintf(&yrot_val[13], "%4.0f", theta[1]);
  yrot_val[17] = ' ';
  sprintf(&zrot_val[13], "%4.0f", theta[2]);
  zrot_val[17] = ' ';
  
  glMatrixMode(GL_PROJECTION);
  glPushMatrix();
    glLoadIdentity();
    glOrtho(0, windoww, 0, windowh, -1, 1);
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
      glLoadIdentity();
      glColor3f(0.0, 0.0, 0.0);
      /* upper left */
      glRasterPos2i(5, windowh-40);
      printstring(GLUT_BITMAP_HELVETICA_18, "pde3_gl");
      glRasterPos2i(100, windowh-35);
      printstring(GLUT_BITMAP_HELVETICA_12, "black is computed, orange is solution");
      glRasterPos2i(5, windowh-60);
      printstring(GLUT_BITMAP_HELVETICA_12, iter_val);
      glRasterPos2i(5, windowh-80);
      printstring(GLUT_BITMAP_HELVETICA_12, err_val);
      glRasterPos2i(5, windowh-100);
      printstring(GLUT_BITMAP_HELVETICA_12, xrot_val);
      glRasterPos2i(5, windowh-120);
      printstring(GLUT_BITMAP_HELVETICA_12, yrot_val);
      glRasterPos2i(5, windowh-140);
      printstring(GLUT_BITMAP_HELVETICA_12, zrot_val);
      glRasterPos2i(5, windowh-160);
      printstring(GLUT_BITMAP_HELVETICA_12, xmin_val);
      glRasterPos2i(5, windowh-180);
      printstring(GLUT_BITMAP_HELVETICA_12, xmax_val);
      glRasterPos2i(5, windowh-200);
      printstring(GLUT_BITMAP_HELVETICA_12, ymin_val);
      glRasterPos2i(5, windowh-220);
      printstring(GLUT_BITMAP_HELVETICA_12, ymax_val);
      glRasterPos2i(5, windowh-240);
      printstring(GLUT_BITMAP_HELVETICA_12, zmin_val);
      glRasterPos2i(5, windowh-260);
      printstring(GLUT_BITMAP_HELVETICA_12, zmax_val);
      /* lower left */
      glRasterPos2i(5, 41);
      printstring(GLUT_BITMAP_HELVETICA_12, "f  or  s changes fast, slow rotation");
      glRasterPos2i(5, 29);
      printstring(GLUT_BITMAP_HELVETICA_12, "W or w changes to wire frame");
      glRasterPos2i(5, 17);
      printstring(GLUT_BITMAP_HELVETICA_12, "arrow keys change rotation axis");
      glRasterPos2i(5, 5);
      printstring(GLUT_BITMAP_HELVETICA_12, "mouse button selects rotation axis");
      glMatrixMode(GL_PROJECTION);
      /* lower right */
      glColor4f(0.5, 1.0, 0.5, 0.0);
      button(windoww-400, 5, windoww-305, 30, "pause");
      glColor4f(0.7, 0.7, 0.7, 1.0);
      button(windoww-300, 5, windoww-205, 30, "Step");
      glColor3f(0.5,0.5,1.0);
      button(windoww-200, 5, windoww-105, 30, "Run");
      glColor3f(1.0,0.5,0.5);
      button(windoww-100, 5, windoww-  5, 30, "Restart");
    glPopMatrix();
    glMatrixMode(GL_MODELVIEW);
  glPopMatrix();
} /* end printhelp */

static void display(void)
{
  /* display callback, clear frame buffer and z buffer, */
  /* rotate _pde and draw, swap buffers */
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glLoadIdentity();
  glRotatef(theta[0], 1.0, 0.0, 0.0);
  glRotatef(theta[1], 0.0, 1.0, 0.0);
  glRotatef(theta[2], 0.0, 0.0, 1.0);
  draw_pde();
  printhelp();
  glFlush();
  glutSwapBuffers();
} /* end display */

static void spin_pde(void)
{
  /* Idle callback, spin _pde 'rotation' degrees about selected axis */
  theta[axis] = theta[axis] + rotation;
  if(theta[axis] > 360.0) theta[axis] = theta[axis] - 360.0;
  if(diff>0.0001 && iter<max_iter && (run || step))
  {
    iterate();
    check();
    avgerror = error/((nx-2)*(ny-2));
    printf("iter iter=%d, diff=%g, error=%g, avg=%g, max=%g\n",
           iter, diff, error, avgerror, maxerror);
    iter++;
    step = 0;
  }
  glutPostRedisplay();
} /* end spin_pde */

static int in_rect(int x, int y, int w, int h)
{
  return xm>x && xm<x+w && ym>y && ym<y+h;
} /* end in_rect */

static void mouse(int btn, int state, int x, int y)
{
  if(btn==GLUT_RIGHT_BUTTON && state == GLUT_DOWN)
  {
    axis++;
    if(axis>2) axis=0;
    if(theta[axis]==0.0) theta[axis]= 0.25;
  }
  xm = x; /* global mouse coordinates */
  ym = windowh-y;
  if(btn==GLUT_LEFT_BUTTON && state==GLUT_DOWN)
  {
    if(in_rect(windoww-100, 5, windoww-5, 30)) /* restart */
    {
      iter = 0;
      init_boundary();
      diff = 1.0e6;
    }
    else if(in_rect(windoww-200, 5, windoww-105, 30)) /* run */
    {
      run = 1;
    }
    else if(in_rect(windoww-300, 5, windoww-205, 30)) /* step */
    {
      step = 1;
    }
    else if(in_rect(windoww-400, 5, windoww-305, 30)) /* pause */
    {
      run = 0;
    }
  }
  glutPostRedisplay();  
} /* end mouse */

static void keyboard(unsigned char key, int x, int y)
{
  if(key == 'f') {rotation=rotation+0.25; if(rotation>5.0) rotation=5.0;}
  if(key == 's') {rotation=rotation-0.25; if(rotation<0.0) rotation=0.0;}
  if(key == 'W') wireframe=1;
  if(key == 'w') wireframe=0;
  if(key == 'p') run = 0;
  if(key == 'r') run = 1;
  /* Use x, X, y, Y, z, Z to change angle about axis */
  if(key == 'x') theta[0] -= 2.0;
  if(key == 'X') theta[0] += 2.0;
  if(key == 'y') theta[1] -= 2.0;
  if(key == 'Y') theta[1] += 2.0;
  if(key == 'z') theta[2] -= 2.0;
  if(key == 'Z') theta[2] += 2.0;
  glutPostRedisplay();
} /* end keyboard */

static void special(int k, int x, int y)
{
  switch(k)
  {
    case GLUT_KEY_LEFT:
      axis = 1;
      break;
    case GLUT_KEY_RIGHT:
      axis = 2;
      break;
    case GLUT_KEY_DOWN:
      axis = 0;
      break;
    case GLUT_KEY_UP:
      wireframe = 1-wireframe;
      break;
  }
} /* end sprcial */

static void myReshape(int w, int h)
{
  glViewport(0, 0, w, h);
  windoww = w;
  windowh = h;
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  if(w <= h)
      glOrtho(-1.5, 1.5, -1.5 * (GLfloat)h / (GLfloat)w,
              1.5 * (GLfloat)h / (GLfloat)w, -10.0, 10.0);
  else
      glOrtho(-1.5 * (GLfloat)w / (GLfloat)h,
              1.5 * (GLfloat)w / (GLfloat)h, -1.5, 1.5, -10.0, 10.0);
  glMatrixMode(GL_MODELVIEW);
  glClearColor(1.0, 1.0, 1.0, 1.0);
} /* end myReshape */

int main(int argc, char *argv[])
{
  printf("pde3.c running\n");
  h = (xmax-xmin)/(nx-1); /* step size in x */
  /*= (ymax-ymin)/(ny-1);    step size in y */
  printf("xmin=%g, xmax=%g, nx=%d, h=%g \n", xmin, xmax, nx, h);
  printf("ymin=%g, ymax=%g, ny=%d, h=%g \n", ymin, ymax, ny, h);
  init_boundary();
  scalex = 2.0/(xmax-xmin);
  scaley = 2.0/(ymax-ymin);
  scalez = 2.0/(zmax-zmin);
  offsetx = -(xmax+xmin)/2.0;
  offsety = -(ymax+ymin)/2.0;
  offsetz = -(zmax+zmin)/2.0;

  glutInit(&argc, argv);
  /* need both double buffering and z buffer */
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
  glutInitWindowSize(windoww, windowh);
  glutCreateWindow(argv[0]);
  glutReshapeFunc(myReshape);      /* enable resize/reshape */
  glutDisplayFunc(display);        /* enable display */
  glutIdleFunc(spin_pde);          /* enable idle function */
  glutMouseFunc(mouse);            /* enable mouse */
  glutKeyboardFunc(keyboard);      /* enable keyboard */
  glutSpecialFunc(special);        /* enable arrow keys */
  glEnable(GL_DEPTH_TEST);         /* enable hidden surface removal */
  glutMainLoop();
  return 0;
} /* end pde3_gl.c */