/* do_tri_split.c  split xx.coord    xx.tri    xx.bound  into
 *                       xxs.coord   xxs.tri   xxs.bound
 *
 * input geometry, three files,  root= xx  .coord,  .tri,  .bound
 * split every triangle x0,y0 x1,y1 x2,y2 into four triangles
 * x01=(x0+x1)/2  y01=(y0+y1)/2  midpoints, also x02 ...  x12 ... 
 * output x0,y0 x01,y01 x02,y02
 * output x1,y1 x01,y01 x12,y12
 * output x2,y2 x02,y02 x12,y12
 * output x01,y01 x02,y02 x12,y12
 */

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>

static int do_input(char fname[]);
static int do_output(char outname[]);
static void uniqueint(int *n, int a[]);
static void sortint(int n, int a[]);
static void freevert(int *nc, int c[], int na, int a[], int nb, int b[]);
static int not_in(int i, int a[], int n);
static void tri_split(void);
static void split_bound(int i1, int i2, int inew);
static int find_coord(double x, double y);

#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, vertices */
#define sz 3600
#define sz3 10800

static int    debug=3;         /* debug print level */
static char   root[60]="22_t"; /* root name for .tri, .coord, .bound */
static char   oroot[60]="22_ts"; /* root name of split files */ 
static int    minvert;         /* minimum vertex number, reduced to zero */
static int    ntri;            /* number of triangles */ 
static int    t1[sz], t2[sz], t3[sz]; /* vertex numbers on triangles */
static int    ncoord;           /* number of coordinates, also equals nvert */
static double xc[sz], yc[sz];   /* coordinates in vertex order */
static int    b1[sz], b2[sz];   /* dirichlet boundary edge vertices */
static int uniquev[sz3];        /* all vertices */
static int nvert;               /* number of unique vertices */
static int uniqueb[sz3];        /* bound vertices */
static int nbound;              /* number of boundary vertices */
static int uniquef[sz3];        /* free vertices, 0 to nfree-1 */
static int nfree;               /* degrees of freedom  ntri-nbound */
static int nuniquev, nuniqueb, nuniquef; /* vertices counts */
static double xmin, xmax, ymin, ymax;

int main(int argc, char *argv[])
{
  int stat;
  char fname[60], oname[60];

  printf("do_split_tri.c running \n");
  printf("Given xx.coord  xx.tri  xx.bound \n");
  printf("Split each triangle into four \n");
  printf("Write xxs.coord  xxs.tri  xxs.bound \n");

  if(argc>2)
  {
    strcpy(fname,argv[1]);
    strcpy(oname,argv[2]);
  }
  else
  {
    strcpy(fname,root);
    strcpy(oname,oroot);
  }
  printf("reading files %s.coord, %s.tri, %s.bound, writing %s \n",
	 fname, fname, fname, oname);

  stat = do_input(fname);
  if(stat != 0) return stat; /* can not continue */

  /* splitting input triangles */
  tri_split();

  strcat(fname,"s");
  stat = do_output(oname);
  if(stat != 0) return stat; /* can not continue */
  printf("triangles=%d, vertices=%d, boundary=%d, free=%d \n",
         ntri, nvert, nuniqueb, nfree); 
  printf("do_tri_split.c finished. \n");
  return 0;
} /* end main of do_tri_split */

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

  /* 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 */
  nuniquev=0;
  for(i=0; i<ntri; i++)
  {
    t1[i]=t1[i]-minvert;
    t2[i]=t2[i]-minvert;
    t3[i]=t3[i]-minvert;
    uniquev[nuniquev]=t1[i];
    uniquev[nuniquev+1]=t2[i];
    uniquev[nuniquev+2]=t3[i];
    nuniquev=nuniquev+3;
  }
  uniqueint(&nuniquev, uniquev);

  /* read in dirichlet boundary, set nbound and nuniqueb */
  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);
  nuniqueb=0;
  for(i=0; i<nbound; i++)
  {
    b1[i]=b1[i]-minvert;
    b2[i]=b2[i]-minvert;
    uniqueb[nuniqueb]=b1[i];
    uniqueb[nuniqueb+1]=b2[i];
    nuniqueb=nuniqueb+2;
  }
  uniqueint(&nuniqueb,uniqueb);
  freevert(&nuniquef, uniquef, nuniquev, uniquev, nuniqueb, uniqueb);
  if(debug>0) for(i=0; i<nuniqueb; i++) printf("unique boundary %d \n", uniqueb[i]);
  if(debug>0) printf("number of free vertices is %d \n", nuniquef);
  if(debug>0) for(i=0; i<nuniquef; i++) printf("free vertex %d \n", uniquef[i]);

  /* 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);
  nfree=nvert-nuniqueb;

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

static int do_output(char root[])
{
  int i;
  FILE * inp;
  int stat;
  char fname[60];

  /* write out triangles,  ntri */
  strcpy(fname, root);
  strcat(fname,".tri");
  if(debug>2) printf("about to write triangles to %s \n", fname);
  inp=fopen(fname, "w");
  if(inp==NULL)
  {
    printf("can not open %s for writing \n", fname);
    return 1;
  }
  for(i=0; i<ntri; i++)
  {
    fprintf(inp, "%d %d %d \n", t1[i], t2[i], t3[i]);
  }
  fclose(inp);


  /* write out boundary, nbound  */
  strcpy(fname, root);
  strcat(fname,".bound");
  if(debug>2) printf("about to write boundary to %s \n", fname);
  inp=fopen(fname, "w");
  if(inp==NULL)
  {
    printf("can not open %s for write \n", fname);
    return 1;
  }
  for(i=0; i<nbound; i++)
  {
    fprintf(inp, "%d %d \n", b1[i], b2[i]);
  }
  fclose(inp);

  /* write out coordinates, nvert */
  strcpy(fname, root);
  strcat(fname,".coord");
  if(debug>2) printf("about to write coordinates to %s \n", fname);
  inp=fopen(fname, "w");
  if(inp==NULL)
  {
    printf("can not open %s for writing \n", fname);
    return 1;
  }
  for(i=0; i<nvert; i++)
  {
    fprintf(inp, "%f %f \n", xc[i], yc[i]);
  }
  fclose(inp);
  return 0;
} /* end do_output */


static void uniqueint(int * n, int a[])
{
  int nold=*n;
  int i, j;
  if(debug>2) printf("finding unique of nold=%d \n", nold);
  sortint(nold,a);
  j=0;
  for(i=1; i<nold; i++)
  {
    if(a[i]>a[i-1])
    {
      j++;
      a[j]=a[i];
    }
  }
  j++;
  *n=j;
  if(debug>2) printf("found unique of n=%d \n", j);
} /* end uniqueint */

static void sortint(int n, int a[])
{
  int i, j, temp;
  for(i=0; i<n-1; i++)
    for(j=i+1; j<n; j++)
      if(a[i]>a[j]) {temp=a[i]; a[i]=a[j]; a[j]=temp;}
} /* end sortint */

static void freevert(int *nc, int c[], int na, int a[], int nb, int b[])
{
  /* a and b sorted integer arrays, c will be from a, not in b */
  int ia=0;
  int ib=0;
  int ic=0;
  if(debug>2) printf("freevert na=%d, nb=%d \n", na, nb);
  while(1)
  {
    if(ib>=nb || a[ia]!=b[ib]) {c[ic]=a[ia]; ic++; ia++;}
    else             {ia++; ib++;}
    if(ia>=na) break;
  }
  if(debug>2) printf("freevert nc free = %d \n", ic);
  *nc=ic;
} /* end freevert */


static int not_in(int i, int a[], int n)
{
  int found=0;
  int j;

  for(j=0; j<n; j++)
  {
    if(i==a[j]) { found=1; break; }
  }
  return 1-found;
} /* end not_in */



/* tri_split  split triangles, one becomes four  */
void tri_split()
{
  int i, j12, j23, j31;
  int nvert_old, ntri_old;
  int i1, i2, i3, i12, i23, i31; /* indices of new split points */
  double x1, y1, x2, y2, x3, y3, x12, y12, x23, y23, x31, y31;

  nvert_old=nvert;
  ntri_old=ntri;
  for(i=0; i<ntri_old; i++)
  {
    i1=t1[i];
    i2=t2[i];
    i3=t3[i];
    x1=xc[i1];
    y1=yc[i1];
    x2=xc[i2];
    y2=yc[i2];
    x3=xc[i3];
    y3=yc[i3];

    j12=find_coord((x1+x2)/2.0,(y1+y2)/2.0); /* may update nvert,xc,yc */
    j23=find_coord((x2+x3)/2.0,(y2+y3)/2.0); /* may update nvert,xc,yc */
    j31=find_coord((x3+x1)/2.0,(y3+y1)/2.0); /* may update nvert,xc,yc */

    t1[i]=i1; t2[i]=j12; t3[i]=j31;
    t1[ntri]=j12; t2[ntri]=i2; t3[ntri]=j23;
    ntri++;
    t1[ntri]=j31; t2[ntri]=j12; t3[ntri]=j23;
    ntri++;
    t1[ntri]=j31; t2[ntri]=j23; t3[ntri]=i3;
    ntri++;
    split_bound(i1, i2, j12);
    split_bound(i2, i3, j23);
    split_bound(i3, i1, j31);
  }
  nuniquev=nvert;
  nfree=nvert-nuniqueb;
  printf("tri_split nvert_old=%d, nvert=%d, ntri_old=%d, ntri=%d \n",
         nvert_old, nvert, ntri_old, ntri);
} /* end tri_split */

static int find_coord(double x, double y)
{
  int i, j;

  for(i=0; i<nvert; i++)
  {
    if(abs(xc[i]-x)+abs(yc[i]-y)<0.0001)
      return i;
  }
  j=nvert;
  xc[j]=x;
  yc[j]=y;
  nvert++;
  return j;
} /* end find coord */

static void split_bound(int i1, int i2, int inew)
{
  int i, jbound;

  /* if i1,i2 boundary, change i1,i2  i1,inew, add inew,i2  */
  jbound = nbound;
  for(i=0; i<jbound; i++)
  {
    if((b1[i]==i1 && b2[i]==i2) || (b1[i]==i2 && b2[i]==i1))
    {
      b2[nbound] = b2[i];
      b2[i] = inew;
      b1[nbound] = inew;
      nbound++;
      break;
    }  
  }

  /*  fix nfree, nuniqueb, uniquev */
  if(not_in(i1, uniqueb, nuniqueb)==0 && not_in(i2, uniqueb, nuniqueb)==0)
  {
    uniqueb[nuniqueb]=inew;
    nuniqueb++;
  }
} /* end split_bound */

/* end do_tri_split.c */