/* 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>
#include "point_in_poly.h" /* only for pitri() point in triangle */

/* 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 nodes4     component names      */
/* noden1  noden2  noden3 noden4     component lengths, sum is n_ndata  */
/* d4 struct {double dirichlet; double fourth; double fifth; double sixth;} */

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

/* model_data is optional fourth dimension for vertices             */
/* modelcomp  models1  models2  models3 models4  component names    */
/* modeln1  modeln2  modeln3  modeln4 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", &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, 2, 3, 4  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;
        noden4 = 0;
        stat = fgets(line, 1023, inp); /* string */
        len = strlen(line);
        nodes1 = (char *)malloc(len);
        strncpy(nodes1, line, len);
        nodes2 = (char *)NULL;
        nodes3 = (char *)NULL;
        nodes4 = (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;
        noden4 = 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;
        nodes4 = (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 if(n_ndata==3)
    {
      stat = fgets(line, 1023, inp);
      sscanf(line,"%d %d %d %d", &nodecomp, &noden1, &noden2, &noden3);
      if(nodecomp==3 && noden1==1 && noden2==1 && noden3==1)
      {
        noden4 = 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);
        stat = fgets(line, 1023, inp); /* string */
        len = strlen(line);
        nodes3 = (char *)malloc(len);
        strncpy(nodes3, line, len);
        nodes4 = (char *)NULL;

        d4th = (d4 *)malloc(sizeof(d4)*(n_vertices+1));
        for(i=1; i<=n_vertices; i++)
        {
          fscanf(inp,"%d %lf %lf %lf\n", &k, &d4th[i].dirichlet,
                 &d4th[i].fourth, &d4th[i].fifth);
        }
      }
      else
      {
        printf("n_data=3, requires 3 1 1 1 for this application, exiting.\n");
        exit(1);
      }
    }
    else if(n_ndata==4)
    {
      stat = fgets(line, 1023, inp);
      sscanf(line,"%d %d %d %d %d", &nodecomp, &noden1, &noden2,
             &noden3, &noden4);
      if(nodecomp==4 && noden1==1 && noden2==1 && noden3==1 && noden4==1)
      {
        noden4 = 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);
        stat = fgets(line, 1023, inp); /* string */
        len = strlen(line);
        nodes3 = (char *)malloc(len);
        strncpy(nodes3, line, len);
        stat = fgets(line, 1023, inp); /* string */
        len = strlen(line);
        nodes4 = (char *)malloc(len);
        strncpy(nodes4, line, len);

        d4th = (d4 *)malloc(sizeof(d4)*(n_vertices+1));
        for(i=1; i<=n_vertices; i++)
        {
          fscanf(inp,"%d %lf %lf %lf %lf\n", &k, &d4th[i].dirichlet,
                 &d4th[i].fourth, &d4th[i].fifth, &d4th[i].sixth);
        }
      }
      else
      {
        printf("n_data=4, requires 4 1 1 1 1 for this application, exiting.\n");
        exit(1);
      }
    }


    else
    {
      printf("n_data not 1, 2, 3, 4 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 */

int pde_closed_ucd() /* called after pde_read_ucd, by user */
{
  int i, j, k, m;
  int unique = 1;
  double tol = 0.001;
  int closed = 1;
  int edge1[10000];
  int edge2[10000];
  int cnt  [10000];
  int nedge = 0;
  int e1, e2, et, new;
  double dif;
  int  ctypn[9] = { 0,      1,     2,       3,     4,
                    4,      5,     6,       8};

  /* check for duplicate vertices  tol = 0.001 */
  for(i=1; i<=n_vertices-1; i++)
  { 
    for(j=i+1; j<=n_vertices; j++)
    { 
      dif = abs(vert[i].x-vert[j].x) +
            abs(vert[i].y-vert[j].y) +
	    abs(vert[i].z-vert[j].z);
      if(dif > tol) unique = 0;
    }
  }

  /* check for every edge in exactly two polygons */
  for(i=1; i<=n_cells; i++)
  {
    k = ctypn[cell[i].ctyp];
    e1 = cell[i].v[k]; // last
    for(j=1; j<=k; j++)
    {
      e2 = cell[i].v[j];
      if(e1>e2) {et = e2; e2 = e1; e1 = et;}
      new = 1;
      for(m=0; m<nedge; m++)
      {
	if(e1==edge1[m] && e2==edge2[m])
	{
	  cnt[m]++;
	  new = 0;
	}
      } // end m
      if(new==1)
      {
        edge1[nedge] = e1;
        edge2[nedge] = e2;
        cnt[nedge] = 1;
        nedge++;
      }
      e1 =  cell[i].v[j];
    } // end edges j in this cell
  } // end i cell
  printf("checking %d edges for cnt==2 \n", nedge);
  for(m=0; m<nedge; m++)
  {
    if(cnt[m]!=2)
    {
      printf("bad edge %d, %d, cnt=%d \n",
	     edge1[m], edge2[m], cnt[m]);
      if(cnt[m]==1) closed = 0;
    }
  }
  if(closed == 0) printf("polyhedra not closed \n");

  return closed + unique; // should be 3
} /* end pde_closed_ucd */

int pde_inside3_ucd(double x, double y, double z)
{
  int c = 0;
  int i, j, k, ct;
  float vertx[8], verty[8], vertz[8];
  int  ctypn[9] = {0, 1, 2, 3, 4, 4, 5, 6, 8};

  for(i=1; i<=n_cells; i++)
  {
    k = ctypn[cell[i].ctyp]; // only valid for k==3
    if(k!=3) printf("non triangle in pde_inside3_ucd i=%d \n", i);
    for(j=1; j<=k; j++)
    {
      vertx[j-1] = vert[cell[i].v[j]].x;
      verty[j-1] = vert[cell[i].v[j]].y;
      vertz[j-1] = vert[cell[i].v[j]].z;
    }
    ct = pitri(3, vertx, verty, vertz, x, y, z);
    c = c + ct;
  }
  return c;
} /* end pde_inside3_ucd */

/* end pde_read_ucd.c */