/* pde_read_ucd.c  all variables and arrays named in pde_read_ucd.h */
/* from  pde_read_ucd.c  just call    pde_read_ucd(file_name);      */
/*                                    only one call per execution   */
/* start with subscript 1, <=n_vertices   etc.                      */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
/* no _MAIN_ GLOBAL==extern                                         */
#include "pde_read_ucd.h"
/* n_vertices  n_cells  n_ndata  n_cdata  n_mdata                   */
/* nodes called vertices in this implementation                     */

/* vertex data  x,y,z  4th optional vertex under model data         */
/* vert struct {double x; double y; double z;}                      */
/* key  struct {int u; int i;}                                      */

/* cell data                                                        */
/* cell struct {int i; int matl; int ctyp; int v[8];}               */
/* ctyp  pt=1(1)   line=2(2) tri=3(3)    quad=4(4)                  */
/*       tet=5(4)  pyr=6(5)  prism=7(6)  hex=8(8)  type number, vert*/

/* node data is Dirichlet boundary values                           */
/* nodecomp  nodes1  nodes2  nodes3     component names             */
/* noden1  noden2  noden3      component lengths, sum is n_ndata    */
/* d4 struct {double dirichlet; double fourth;}                     */

/* cell_data is Neumann boundary values, outward normal positive    */
/* cellcomp  cells1, cells2, cells3     component names             */
/* celln1  celln2  celln3       component lengths, sum is n_cdata   */

/* model_data is optional fourth dimension for vertices             */
/* modelcomp  models1  models2  models3   component names           */
/* modeln1  modeln2  modeln3 component lengths, sum is n_mdata      */
/*                     for this application      n_mdata=n_vertices */

void pde_read_ucd(char * file_name)
{
  int i, j, k, len;
  FILE * inp;
  char line[1024];
  char *stat = NULL;
  float x,y,z;
  int icell, matl;
  char typ[8];
  /*  cell type name, integer number, number of vertices */
  int cvert[8];
  char *ctyps[9] = {"none", "pt",  "line",  "tri", "quad",
                   "tet",  "pyr", "prism", "hex"};
  int  ctypi[9] = { 0,      1,     2,       3,     4,
                    5,      6,     7,       8};
  int  ctypn[9] = { 0,      1,     2,       3,     4,
                    4,      5,     6,       8};


  inp = fopen(file_name,"r");
  if(inp==NULL)
  {
    printf("file %s could not be opened for reading by pde_read_ucd\n",
           file_name);
    exit(1);
  }
  while(1) /* read comments from UCD file */
  {  
    stat = fgets(line, 1023, inp);
    if(stat==NULL) 
    {
      printf("pde_read_ucd, premature EOF in %s\n", file_name);
      exit(1);
    }
    /* printf("%s",line); */ /* optional */
    if(line[0]!='#') break;
  } /* end while to get past comments */
  sscanf(line,"%d %d %d %d %d", &n_vertices, &n_cells,
         &n_ndata, &n_cdata, &n_mdata);

  /* vertices = nodes */
  vert = (xyz *)malloc(sizeof(xyz)*(n_vertices+1));
  key = (keys *)malloc(sizeof(keys)*(n_vertices+1));
  for(i=1; i<=n_vertices; i++)
  {
    key[i].i = i;
    stat = fgets(line, 1023, inp);
    sscanf(line,"%d %f %f %f\n", &key[i].u, &x, &y, &z);
    vert[i].x=x;
    vert[i].y=y;
    vert[i].z=z;
  } /* end vertices */

  /* cells */
  cell = (cells *)malloc(sizeof(cells)*(n_cells+1));
  for(i=1; i<=n_cells; i++)
  {
    stat = fgets(line, 1023, inp);
    sscanf(line,"%d %d %s", &icell, &matl, typ);
    cell[i].i = icell; /* not changed */  
    cell[i].matl = matl;
    k = 0;
    for(j=1; j<9; j++)
    {
      if(strcmp(ctyps[j],typ)==0)
      {
        cell[i].ctyp = ctypi[j];
        k = ctypn[j];
        break;
      }
    }
    switch(k)
    {
    case 1:
      sscanf(line,"%d %d %s %d", &icell, &matl, typ, &cell[i].v[1]);
      break;
    case 2:
      sscanf(line,"%d %d %s %d %d", &icell, &matl, typ, &cell[i].v[1],
             &cell[i].v[2]);
      break;
    case 3:
      sscanf(line,"%d %d %s %d %d %d", &icell, &matl, typ, &cell[i].v[1],
             &cell[i].v[2], &cell[i].v[3]);
      break;
    case 4:
      sscanf(line,"%d %d %s %d %d %d %d", &icell, &matl, typ, &cell[i].v[1],
             &cell[i].v[2], &cell[i].v[3], &cell[i].v[4]);
      break;
    case 5:
      sscanf(line,"%d %d %s %d %d %d %d %d", &icell, &matl, typ, &cell[i].v[1],
             &cell[i].v[2], &cell[i].v[3], &cell[i].v[4], &cell[i].v[5]);
      break;
    case 6:
      sscanf(line,"%d %d %s %d %d %d %d %d %d", &icell, &matl, typ,
             &cell[i].v[1], &cell[i].v[2], &cell[i].v[3], &cell[i].v[4],
             &cell[i].v[5], &cell[i].v[6]);
      break;
    case 8:
      sscanf(line,"%d %d %s %d %d %d %d %d %d %d %d %d", &icell, &matl, typ,
             &cell[i].v[1], &cell[i].v[2], &cell[i].v[3], &cell[i].v[4],
             &cell[i].v[5], &cell[i].v[6], &cell[i].v[7], &cell[i].v[8]);
      break;
    default:
      printf("bad cell type %s\n", line);
    } /* end switch */
  } /* end cells */

  /* node=vertex data Dirichlet boundary values */
  if(n_ndata>0) /* must be zero, 1, or 2 for this application */
  {
    if(n_ndata==1)
    {
      stat = fgets(line, 1023, inp);
      sscanf(line,"%d %d", &nodecomp, &noden1);
      if(nodecomp==1 && noden1==1)
      {
        noden2 = 0;
        noden3 = 0;
        stat = fgets(line, 1023, inp); /* string */
        len = strlen(line);
        nodes1 = (char *)malloc(len);
        strncpy(nodes1, line, len);
        nodes2 = (char *)NULL;
        nodes3 = (char *)NULL;
        dirichlet_bound = (double *)malloc(sizeof(double)*(n_vertices+1));
        for(i=1; i<=n_vertices; i++)
        {
          fscanf(inp,"%d %lf\n", &k, &dirichlet_bound[i]);
        }
      }
      else
      {
        printf("n_data=1, requires 1 1 for this application, exiting.\n");
        exit(1);
      }
    }
    else if(n_ndata==2)
    {
      stat = fgets(line, 1023, inp);
      sscanf(line,"%d %d %d", &nodecomp, &noden1, &noden2);
      if(nodecomp==2 && noden1==1 && noden2==1)
      {
        noden3 = 0;
        stat = fgets(line, 1023, inp); /* string */
        len = strlen(line);
        nodes1 = (char *)malloc(len);
        strncpy(nodes1, line, len);
        stat = fgets(line, 1023, inp); /* string */
        len = strlen(line);
        nodes2 = (char *)malloc(len);
        strncpy(nodes2, line, len);
        nodes3 = (char *)NULL;
        d4th = (d4 *)malloc(sizeof(d4)*(n_vertices+1));
        for(i=1; i<=n_vertices; i++)
        {
          fscanf(inp,"%d %lf %lf\n", &k, &d4th[i].dirichlet,
                 &d4th[i].fourth);
        }
      }
      else
      {
        printf("n_data=2, requires 2 1 1 for this application, exiting.\n");
        exit(1);
      }
    }
    else
    {
      printf("n_data not 1 or 2, required for this application, exiting.\n");
      exit(1);
    }

  } /* end node data, Dirichlet values */

  /* cell data  Neumann Boundary values */
  if(n_cdata>0) /* must be zero or 1 for this application */
  {
    stat = fgets(line, 1023, inp);
    sscanf(line,"%d %d", &cellcomp, &celln1);
    celln2 = 0;
    celln3 = 0;
    if(n_cdata!=1 || cellcomp!=1 || celln1!=1)
    {
      printf("non 1, n_cdata or cellcomp or celln1, exiting.\n");
      exit(1);
    }
    stat = fgets(line, 1023, inp); /* string */
    len = strlen(line);
    cells1 = (char *)malloc(len);
    strncpy(cells1, line, len);
    cells2 = (char *)NULL;
    cells3 = (char *)NULL;
    neumann_bound = (double *)malloc(sizeof(double)*(n_cells+1));
    for(i=1; i<=n_cells; i++)
    {
      fscanf(inp,"%d %lf\n", &k, &neumann_bound[i]);
    }
  } /* end cell data, Dirichlet values */

  /* model data, optional 4th vertex component for 4D */
  if(n_mdata>0)  /* must be zero or 1 for this application */
  {
    stat = fgets(line, 1023, inp);
    sscanf(line,"%d %d", &modelcomp, &modeln1);
    modeln2 = 0;
    modeln3 = 0;
    if(n_mdata!=1 || modelcomp!=1 || modeln1!=1)
    {
      printf("non 1, n_mdata or modelcomp or modeln1, exiting.\n");
      exit(1);
    }
    stat = fgets(line, 1023, inp); /* string */
    len = strlen(line);
    models1 = (char *)malloc(len);
    strncpy(models1, line, len);
    models2 = (char *)NULL;
    models3 = (char *)NULL;
    fourth_dim = (double *)malloc(sizeof(double)*(n_vertices+1));
    for(i=1; i<=n_vertices; i++)
    { 
      fscanf(inp,"%d %lf\n", &k, &fourth_dim[i]);
    }
  } /* end model data, fourth coordinate for vertices, 4D PDE */
  /* just ignore extra stuff in file */
  fclose(inp);
} /* end pde_read_ucd.c */