/* plot_fem_tri.c  C.coord, C.tri, C.bound */

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

static int do_input(); /* .coord, .tri, .bound */

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

/* maximum number of triangles or vertices, *3, squared */
#define sz 200
#define sz3 600

static int    debug=3;          /* debug print level */
static char   root[]="C";       /* root name for .tri, .coord, .bound */
static int    ntri;             /* number of triangles */ 
static int    minvert;          /* minimum vertex number, reduced to zero */
static int    nvert;            /* number of vertices */
static int    t1[sz], t2[sz], t3[sz]; /* vertex numbers on triangles */
static int    ncoord;           /* number of coordinates, also equals nvert */
static double xc[sz3], yc[sz3]; /* coordinates in vertex order */
static int    nbound;           /* number of boundary pairs */
static int    b1[sz], b2[sz];   /* dirichlet boundary edge vertices */
static double xmin, xmax, ymin, ymax; /* determined by input */

static int width=500;
static int height=500;

void display(void)
{
  char text[]="plot_fem_tri";
  char *p;
  int i, j, k;
  double x1, y1, x2, y2, x3, y3;
  float cred, cgreen, cblue, fct;
  
  /* clear window */
  glClear(GL_COLOR_BUFFER_BIT); 
  glLoadIdentity ();
  glColor3f(0.0, 0.0, 0.0);

  /* draw text, in its own context */
  glPushMatrix();
    glLoadIdentity ();
    glColor3f(0.0, 0.0, 0.0);
    glEnable(GL_LINE_SMOOTH);
    glTranslatef(-0.95, 0.89, 0.0);
    glScalef(0.00055, 0.00055, 0.0);
    for(p=text; *p; p++)
      glutStrokeCharacter(GLUT_STROKE_ROMAN, *p);
  glPopMatrix();
  
  /* Draw the segments */
  glColor3f(0.0, 0.0, 0.0);
  glBegin(GL_LINES);
    for(i=0; i<ntri; i++)
    {
      x1=xc[t1[i]];
      y1=yc[t1[i]];
      x2=xc[t2[i]];
      y2=yc[t2[i]];
      x3=xc[t3[i]];
      y3=yc[t3[i]];
      glVertex3d(x1, y1, 0.0);
      glVertex3d(x2, y2, 0.0);
      glVertex3d(x2, y2, 0.0);
      glVertex3d(x3, y3, 0.0);
      glVertex3d(x3, y3, 0.0);
      glVertex3d(x1, y1, 0.0);
    }
  glEnd();
  /* draw the boundary points */
  glColor3f(1.0, 0.0, 0.0);
  glPointSize(6.0);
  glBegin(GL_POINTS);
    for(i=0; i<nbound; i++)
    {
      x1=xc[b1[i]];
      y1=yc[b1[i]];
      x2=xc[b2[i]];
      y2=yc[b2[i]];
      glVertex3d(x1, y1, 0.0);
      glVertex3d(x2, y2, 0.0);
    }
  glEnd();
  /* draw the boundary lines */
  glColor3f(0.0, 0.0, 1.0);
  glBegin(GL_LINES);
    for(i=0; i<nbound; i++)
    {
      x1=xc[b1[i]];
      y1=yc[b1[i]];
      x2=xc[b2[i]];
      y2=yc[b2[i]];
      glVertex3d(x1, y1, 0.0);
      glVertex3d(x2, y2, 0.0);
    }
  glEnd();
  glFlush(); 
}


/* This routine handles window resizes */
void reshape(int w, int h)
{
  width = w;
  height = h;
  /* Set the transformations */
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glOrtho(xmin, xmax, ymin, ymax, -1.0, 1.0);
  glMatrixMode(GL_MODELVIEW);
  glViewport(0, 0, w, h);
}

void init()
{
  int stat;
  double dx, dy;

  stat=do_input(); /* read .coord, .tri, .bound files */
  glClearColor(1.0, 1.0, 1.0, 0.0);
  dx=(xmax-xmin)/20.0;
  dy=(ymax-ymin)/20.0;
  xmax = xmax+dx;
  xmin = xmin-dx;
  ymax = ymax+dy;
  ymin = ymin-dy;
  glMatrixMode (GL_PROJECTION);
  glLoadIdentity ();
  glOrtho(xmin, xmax, ymin, ymax, -1.0, 1.0);
  glMatrixMode (GL_MODELVIEW);
  glViewport(0, 0, width, height);
}

int main(int argc, char* argv[])
{
  glutInit(&argc,argv);
  glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);  
  glutInitWindowSize(width, height);
  glutInitWindowPosition(0,0); 
  
  glutCreateWindow(argv[0]); 
  glutDisplayFunc(display);
  glutReshapeFunc(reshape);
  init();
  glutMainLoop();
  return 0;
} /* end main of plot_fem_tri.c */

static int do_input()
{
  int i, j, ii;
  char fname[60];
  FILE * inp;
  int stat;

  /* read in triangles, set ntri */
  strcpy(fname, root);
  strcat(fname,".tri");
  if(debug>2) printf("about to read triangles from %s \n", fname);
  inp=fopen(fname, "r");
  if(inp==NULL)
  {
    printf("can not open %s for reading \n", fname);
    return 1;
  }
  stat=0;
  ntri=0;
  minvert=999999;
  nvert=0;
  if(debug>0) printf("triangles read from %s \n",fname);
  while(stat>=0)
  {
    stat=fscanf(inp, "%d %d %d", &t1[ntri], &t2[ntri], &t3[ntri]);
    if(stat<0) break;
    if(debug>0) printf("tri %2d has vertices %2d  %2d  %2d \n",
                ntri, t1[ntri], t2[ntri], t3[ntri]);
    if(t1[ntri]>nvert) nvert=t1[ntri];
    if(t2[ntri]>nvert) nvert=t2[ntri];
    if(t3[ntri]>nvert) nvert=t3[ntri];
    if(t1[ntri]<minvert) minvert=t1[ntri];
    if(t2[ntri]<minvert) minvert=t2[ntri];
    if(t3[ntri]<minvert) minvert=t3[ntri];
    ntri++;
  }
  fclose(inp);
  printf("subtracting minvert=%d from all vertices, using base zero\n",
         minvert);
  nvert=nvert+1-minvert; /* still must be dense sequence of numbers */
  for(i=0; i<ntri; i++)
  {
    t1[i]=t1[i]-minvert;
    t2[i]=t2[i]-minvert;
    t3[i]=t3[i]-minvert;
  }

  /* read in dirichlet boundary, set nbound */
  strcpy(fname, root);
  strcat(fname,".bound");
  if(debug>2) printf("about to read boundary from %s \n", fname);
  inp=fopen(fname, "r");
  if(inp==NULL)
  {
    printf("can not open %s for reading \n", fname);
    return 1;
  }
  stat=0;
  nbound=0;
  if(debug>0) printf("Dirichlet boundaries from %s \n",fname);
  while(stat>=0)
  {
    stat=fscanf(inp, "%d %d", &b1[nbound], &b2[nbound]);
    if(stat<0) break;
    if(debug>0) printf("boundary segment %2d has vertices %2d  %2d \n",
                nbound, b1[nbound], b2[nbound]);
    if(b1[nbound]>nvert) printf("bad boundary vertex %d \n", b1[nbound]);
    if(b2[nbound]>nvert) printf("bad boundary vertex %d \n", b2[nbound]);
    if(b1[nbound]<minvert) printf("bad boundary vertex %d \n", b1[nbound]);
    if(b2[nbound]<minvert) printf("bad boundary vertex %d \n", b2[nbound]);
    nbound++;
  }
  fclose(inp);
  printf("subtracting minvert=%d from all vertices, using base zero\n", minvert);
  for(i=0; i<nbound; i++)
  {
    b1[i]=b1[i]-minvert;
    b2[i]=b2[i]-minvert;
  }

  /* read in coordinates, check consistency */
  strcpy(fname, root);
  strcat(fname,".coord");
  if(debug>2) printf("about to read coordinates from %s \n", fname);
  inp=fopen(fname, "r");
  if(inp==NULL)
  {
    printf("can not open %s for reading \n", fname);
    return 1;
  }
  stat=0;
  ncoord=0;
  if(debug>0) printf("coordinates read from %s \n",fname);
  while(stat>=0)
  {
    stat=fscanf(inp, "%lf %lf", &xc[ncoord], &yc[ncoord]);
    if(stat<0) break;
    if(debug>0) printf("coordinate %2d at  %5.2f , %5.2f \n",
                       ncoord, xc[ncoord], yc[ncoord]);
    if(ncoord==0)
    {
      xmax=xc[ncoord];
      xmin=xc[ncoord];
      ymax=yc[ncoord];
      ymin=yc[ncoord];
    }
    if(xc[ncoord]>xmax) xmax=xc[ncoord];
    if(xc[ncoord]<xmin) xmin=xc[ncoord];
    if(yc[ncoord]>ymax) ymax=yc[ncoord];
    if(yc[ncoord]<ymin) ymin=yc[ncoord];
    ncoord++;
  }
  fclose(inp);
  if(ncoord!=nvert)
     printf("\n**** wrong number of coordinates=%d, should be %d\n",
	    ncoord, nvert);

  printf("xmin=%g, xmax=%g, ymin=%g, ymax=%g \n", xmin, xmax, ymin, ymax);
  printf("nvert=%d, nbound=%d, ntri=%d \n",
	 nvert, nbound, ntri);
  return 0;
} /* end do_input */