/* test_split_tetra.c                              */
/*      8 tetrahedrons built from each tetrahedron */

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

static float vert[8193][3] = {{0.0, 0.0, 0.0}, /* unused, start with [1] */
       {-1.0,-1.0,-1.0}, {-0.4, 0.9, -0.8},    /* 1,2 */
       { 0.7, 0.4,-1.0}, { 0.5, 0.5, 1.0}};    /* 3,4 */
static int nvert = 4; /* increment before inserting, count not including zero*/
static int ivert[8193][4] = {{1,2,3,4}};
static int npoly = 1; /* number of 4 vertex tetrahedrons in ivert*/
static int bound[8193][3] = {{1,2,3}, {1,2,4}, {1,3,4}, {2,3,4}};
static int nbound = 4;
static int uniqueb[8193] = {1, 2, 3, 4};
static int nuniqueb = 4;
static int nv;
static FILE* ofile;
static char filename[] = "33_ts.tetra";

int in_bound(int k)
{
  int i, found = 0;
  for(i=0; i<nuniqueb; i++)
  {
    if(uniqueb[i]==k) {found=1; break;}
  }
  return found; /* 1 means found */
} /* end in_bound */

void add_ubound(int k)
{
  int i;
  for(i=0; i<nuniqueb; i++)
  {
    if(uniqueb[i]==k) return;
  }
  uniqueb[nuniqueb] = k;
  nuniqueb++;
} /* end add_ubound */

void add_bound(int p1, int p2, int p3)
{
  bound[nbound][0] = p1;
  bound[nbound][1] = p2;
  bound[nbound][2] = p3;
  nbound++;
} /* end add_ubound */

void next_tetra()
{
  int tvert[8193][4];
  int tpoly;
  int ipoly, j, k;
  int p1, p2, p3, p4, p12, p13, p14, p23, p24, p34, p1234;
  
  tpoly = npoly;
  for(ipoly=0; ipoly<tpoly; ipoly++) /* move old polygons to temp */
    for(j=0; j<4; j++)
      tvert[ipoly][j] = ivert[ipoly][j];

  npoly = 0;
  for(ipoly=0; ipoly<tpoly; ipoly++) /* generate 8 new for each old */
  {
    p1 = tvert[ipoly][0];
    p2 = tvert[ipoly][1];
    p3 = tvert[ipoly][2];
    p4 = tvert[ipoly][3];
    nvert++;
    p12 = nvert;
    if(in_bound(p1) && in_bound(p2)) add_ubound(p12);
    vert[nvert][0] = (vert[p1][0]+vert[p2][0])/2.0;
    vert[nvert][1] = (vert[p1][1]+vert[p2][1])/2.0;
    vert[nvert][2] = (vert[p1][2]+vert[p2][2])/2.0;
    nvert++;
    p13 = nvert;
    if(in_bound(p1) && in_bound(p3)) add_ubound(p13);
    vert[nvert][0] = (vert[p1][0]+vert[p3][0])/2.0;
    vert[nvert][1] = (vert[p1][1]+vert[p3][1])/2.0;
    vert[nvert][2] = (vert[p1][2]+vert[p3][2])/2.0;
    nvert++;
    p14 = nvert;
    if(in_bound(p1) && in_bound(p4)) add_ubound(p14);
    vert[nvert][0] = (vert[p1][0]+vert[p4][0])/2.0;
    vert[nvert][1] = (vert[p1][1]+vert[p4][1])/2.0;
    vert[nvert][2] = (vert[p1][2]+vert[p4][2])/2.0;
    nvert++;
    p23 = nvert;
    if(in_bound(p2) && in_bound(p3)) add_ubound(p23);
    vert[nvert][0] = (vert[p2][0]+vert[p3][0])/2.0;
    vert[nvert][1] = (vert[p2][1]+vert[p3][1])/2.0;
    vert[nvert][2] = (vert[p2][2]+vert[p3][2])/2.0;
    nvert++;
    p24 = nvert;
    if(in_bound(p2) && in_bound(p4)) add_ubound(p24);
    vert[nvert][0] = (vert[p2][0]+vert[p4][0])/2.0;
    vert[nvert][1] = (vert[p2][1]+vert[p4][1])/2.0;
    vert[nvert][2] = (vert[p2][2]+vert[p4][2])/2.0;
    nvert++;
    p34 = nvert;
    if(in_bound(p3) && in_bound(p4)) add_ubound(p34);
    vert[nvert][0] = (vert[p3][0]+vert[p4][0])/2.0;
    vert[nvert][1] = (vert[p3][1]+vert[p4][1])/2.0;
    vert[nvert][2] = (vert[p3][2]+vert[p4][2])/2.0;
    nvert++;
    p1234 = nvert; /* center */
    vert[nvert][0] = (vert[p1][0]+vert[p2][0]+
                      vert[p3][0]+vert[p4][0])/4.0;
    vert[nvert][1] = (vert[p1][1]+vert[p2][1]+
                      vert[p3][1]+vert[p4][1])/4.0;
    vert[nvert][2] = (vert[p1][2]+vert[p2][2]+
                      vert[p3][2]+vert[p4][2])/4.0;

    nbound = 0; /* all new from split */
    if(in_bound(p1 ) && in_bound(p12) && in_bound(p13)) add_bound(p1,p12,p13);
    if(in_bound(p1 ) && in_bound(p12) && in_bound(p14)) add_bound(p1,p12,p14);
    if(in_bound(p1 ) && in_bound(p13) && in_bound(p14)) add_bound(p1,p13,p14);
    ivert[npoly][0] = p1;
    ivert[npoly][1] = p12;
    ivert[npoly][2] = p13;
    ivert[npoly][3] = p14;
    npoly++;
    if(in_bound(p2 ) && in_bound(p12) && in_bound(p23)) add_bound(p2,p12,p23);
    if(in_bound(p2 ) && in_bound(p12) && in_bound(p24)) add_bound(p2,p12,p24);
    if(in_bound(p2 ) && in_bound(p23) && in_bound(p24)) add_bound(p2,p23,p24);
    ivert[npoly][0] = p2;
    ivert[npoly][1] = p12;
    ivert[npoly][2] = p23;
    ivert[npoly][3] = p24;
    npoly++;
    if(in_bound(p3 ) && in_bound(p13) && in_bound(p23)) add_bound(p3,p13,p23);
    if(in_bound(p3 ) && in_bound(p13) && in_bound(p34)) add_bound(p3,p13,p34);
    if(in_bound(p3 ) && in_bound(p23) && in_bound(p34)) add_bound(p3,p23,p34);
    ivert[npoly][0] = p3;
    ivert[npoly][1] = p13;
    ivert[npoly][2] = p23;
    ivert[npoly][3] = p34;
    npoly++;
    if(in_bound(p4 ) && in_bound(p14) && in_bound(p24)) add_bound(p4,p14,p24);
    if(in_bound(p4 ) && in_bound(p14) && in_bound(p34)) add_bound(p4,p14,p34);
    if(in_bound(p4 ) && in_bound(p24) && in_bound(p34)) add_bound(p4,p24,p34);
    ivert[npoly][0] = p4;
    ivert[npoly][1] = p14;
    ivert[npoly][2] = p24;
    ivert[npoly][3] = p34;
    npoly++;
    ivert[npoly][0] = p1234;
    ivert[npoly][1] = p12;
    ivert[npoly][2] = p13;
    ivert[npoly][3] = p14;
    npoly++;
    ivert[npoly][0] = p1234;
    ivert[npoly][1] = p12;
    ivert[npoly][2] = p23;
    ivert[npoly][3] = p24;
    npoly++;
    ivert[npoly][0] = p1234;
    ivert[npoly][1] = p13;
    ivert[npoly][2] = p23;
    ivert[npoly][3] = p34;
    npoly++;
    ivert[npoly][0] = p1234;
    ivert[npoly][1] = p14;
    ivert[npoly][2] = p24;
    ivert[npoly][3] = p34;
    npoly++;
  }
} /* end next_tetra */


int main(int argc, char *argv[])
{
  int i;

  /* make output files  33_ts.tetra  .coord  .bound */
  printf("test_split_tetra.c running \n");
  printf("initial tetrahedron\n");
  printf("npoly=%d, nvert=%d, nbound=%d, nuniqueb=%d\n",
         npoly, nvert, nbound, nuniqueb);
  printf("tetrahedrons\n");
  for(i=0; i<npoly; i++)
  {
    printf("% 4d %4d %4d %4d\n",
           ivert[i][0], ivert[i][1], ivert[i][2], ivert[i][3]);
  }
  printf("coordinates\n");
  for(i=1; i<=nvert; i++) /* [0] unused */
  {
    printf("%5i  %12.9f %12.9f %12.9f\n",
           i, vert[i][0], vert[i][1], vert[i][2]);
  }
  printf("boundary faces\n");
  for(i=0; i<nbound; i++)
  {
    printf(" %4d %4d %4d\n",
           bound[i][0], bound[i][1], bound[i][2]);
  }
  printf("boundary vertices\n");
  for(i=0; i<nuniqueb; i++)
  {
    printf(" %4d\n", uniqueb[i]);
  }

  next_tetra();
  printf("\n first split\n");
  printf("npoly=%d, nvert=%d, nbound=%d, nuniqueb=%d\n",
         npoly, nvert, nbound, nuniqueb);
  printf("tetrahedrons\n");
  ofile = fopen(filename,"w");
  for(i=0; i<npoly; i++)
  {
    fprintf(ofile, "%4d %4d %4d %4d\n",
            ivert[i][0], ivert[i][1], ivert[i][2], ivert[i][3]);
    printf("% 4d %4d %4d %4d\n",
           ivert[i][0], ivert[i][1], ivert[i][2], ivert[i][3]);
  }
  fclose(ofile);
  filename[6]='c'; filename[7]='o'; filename[8]='o';
  filename[9]='r'; filename[10]='d'; 
  printf("coordinates\n");
  ofile = fopen(filename,"w");
  for(i=1; i<=nvert; i++) /* [0] unused */
  {
    fprintf(ofile, "%12.9f %12.9f %12.9f\n",
            vert[i][0], vert[i][1], vert[i][2]);
    printf("%5i  %12.9f %12.9f %12.9f\n",
           i, vert[i][0], vert[i][1], vert[i][2]);
  }
  fclose(ofile);
  filename[6]='b'; filename[7]='o'; filename[8]='u';
  filename[9]='n'; filename[10]='d'; 
  ofile = fopen(filename,"w");
  printf("boundary faces\n");
  for(i=0; i<nbound; i++)
  {
    fprintf(ofile, "%4d %4d %4d\n",
            bound[i][0], bound[i][1], bound[i][2]);
    printf(" %4d %4d %4d\n",
           bound[i][0], bound[i][1], bound[i][2]);
  }
  fclose(ofile);
  filename[6]='u'; filename[7]='b'; filename[8]='o';
  filename[9]='u'; filename[10]='n'; 
  ofile = fopen(filename,"w");
  printf("boundary vertices\n");
  for(i=0; i<nuniqueb; i++)
  {
    fprintf(ofile, "%4d\n", uniqueb[i]);
    printf(" %4d\n", uniqueb[i]);
  }
  fclose(ofile);

  printf("\n next split\n");
  filename[4]='t';
  filename[6]='t'; filename[7]='e'; filename[8]='t';
  filename[9]='r'; filename[10]='a'; 
  next_tetra();
  printf("npoly=%d, nvert=%d, nbound=%d, nuniqueb=%d\n",
         npoly, nvert, nbound, nuniqueb);
  printf("tetrahedrons\n");
  ofile = fopen(filename,"w");
  for(i=0; i<npoly; i++)
  {
    fprintf(ofile, "%4d %4d %4d %4d\n",
            ivert[i][0], ivert[i][1], ivert[i][2], ivert[i][3]);
    printf(" %4d %4d %4d %4d\n",
           ivert[i][0], ivert[i][1], ivert[i][2], ivert[i][3]);
  }
  fclose(ofile);
  filename[6]='c'; filename[7]='o'; filename[8]='o';
  filename[9]='r'; filename[10]='d'; 
  printf("coordinates\n");
  ofile = fopen(filename,"w");
  for(i=1; i<=nvert; i++) /* [0] unused */
  {
    fprintf(ofile, "%12.9f %12.9f %12.9f\n",
            vert[i][0], vert[i][1], vert[i][2]);
    printf("%5i  %12.9f %12.9f %12.9f\n",
           i, vert[i][0], vert[i][1], vert[i][2]);
  }
  fclose(ofile);
  filename[6]='b'; filename[7]='o'; filename[8]='u';
  filename[9]='n'; filename[10]='d'; 
  printf("boundary faces\n");
  ofile = fopen(filename,"w");
  for(i=0; i<nbound; i++)
  {
    fprintf(ofile, "%4d %4d %4d\n",
            bound[i][0], bound[i][1], bound[i][2]);
    printf(" %4d %4d %4d\n",
           bound[i][0], bound[i][1], bound[i][2]);
  }
  fclose(ofile);
  filename[6]='u'; filename[7]='b'; filename[8]='o';
  filename[9]='u'; filename[10]='n'; 
  printf("boundary vertices\n");
  ofile = fopen(filename,"w");
  for(i=0; i<nuniqueb; i++)
  {
    fprintf(ofile, "%4d\n", uniqueb[i]);
    printf(" %4d\n", uniqueb[i]);
  }
fclose(ofile);
  return 0;
} /* end test_split_tetra.c */