// tetra_split1.java           split simple tetra into 12 tetrahedrons
//                             11 vertex 48 triangles
//    reads  .dat with 4 points, 4 triangular surfaces
//    writes .dat
//
// needs Matrix.determinant
// needs datread.java compiled
// needs datwrite.java compiled
//
//                 4                    input: 4 triangles
//                /\                    1 2 4
//               / |\     * 11 center   1 2 4
//              /  | \                  2 3 4
//          7  *   |  * 10 midpoint     3 1 4
//            / 6  |   \   
//           /__*__|__ _\               output 8 tetrahedrons
//        1  \     |    / 3             1  5  6  7 not drawn.
//            \    *9  /                2  8  9  5 Each written
//             \   |  /                 3 10  6  8 as 4 triangles
//            5 *  | * 8                4  7  9 10 in  .dat  file.
//               \ |/                   5  6  7 11
//                \/                    8  9  5 11
//                 2                   10  6  8 11
//                                      7  9 10 11
//                                      5  7  9 11
//                                      8  9 10 11
//                                      6  7 10 11
//                                      5  6  8 11
//

public class tetra_split1
{
  int vert[][] = new int[65][3]; //  triangle faces, 4 per tetra,index 1 offset
  int nvert = 0;
  int mnvert = 4;
  int ncntr = 8;
  double tv[][]   = new double[22][3]; // x,y,z for vertices, 1 offset
                                       // 4 midpoints, center
  int ntv = 0;
  double points[][] = new double[65][3]; // x,y,z 0 offset
  double volume[] = new double[22];
  double volumet = 0;
  double volumetot = 0;
  double area = 0.0;
  double vol4;
  int polys[] = new int[400]; // at least 4*vert
  int kp = 0;
  int jj = 0;
  int tet[][] = new int[49][4]; // 4 vert for each tetra
  int ntet = 0;
  datread DR; // = new datread(filename); // defines data structures
  int num_points = 0; // set from datread, reset for datwrite
  int num_polys = 0; // set from datread, reset for datwrite
  int new_points = 0; // 3 added for each poly (triangle)
  int new_polys = 0; // 1 changed and 3 added for each poly
  double xyz[][] = new double[22][3]; // x, y, z for output
  int newpts[] = new int[10000]; // 3, xyz1, xyz2, xyz3 for output
	
  public tetra_split1(String filein, String fileout) // .dat  .dat
  {
    int kk, pt, pts;
    
    System.out.println("tetra_split1.java running");
    System.out.println("datread.java reading "+filein);
    DR = new datread(filein); // defines data structures
    num_points = DR.num_points;
    num_polys = DR.num_polys;
    System.out.println("num_points="+DR.num_points+
			   ", num_polys="+DR.num_polys);
    System.out.println("print vertices");
    for(int i=0; i<num_points; i++) // vertices
    {
      System.out.println(+(i+1)+"  "+DR.data_points[i].x+
			 "  "+DR.data_points[i].y+"  "+DR.data_points[i].z);
    } // end print vertices

    System.out.println(" ");
    System.out.println("print polygons");
    kk = 0;
    for(int i=0; i<num_polys; i++)
    {
      pts = DR.data_polys[kk++];
      System.out.print(+pts); // no cr
      for(int j=0; j<pts; j++)
      {
	pt = DR.data_polys[kk++];
	System.out.print(" "+ pt);
      } // end j
      System.out.println(" ");
    } // end i
    System.out.println(" ");

    if(num_points != 4 || num_polys != 4)
    {
      System.out.println("must have exactly 4 vertices,");
      System.out.println("must have exactly 4 faces, bye");
      System.exit(0);	
    }
    System.out.println(" ");
    tv[0][0] = 0.0; // unused
    tv[0][1] = 0.0;
    tv[0][2] = 0.0;
    // coordinates of vertex
    tv[1][0] = DR.data_points[0].x; // vert 1
    tv[1][1] = DR.data_points[0].y;
    tv[1][2] = DR.data_points[0].z;
    tv[2][0] = DR.data_points[1].x; // vert 2
    tv[2][1] = DR.data_points[1].y;
    tv[2][2] = DR.data_points[1].z;
    tv[3][0] = DR.data_points[2].x; // vert 3
    tv[3][1] = DR.data_points[2].y;
    tv[3][2] = DR.data_points[2].z;
    tv[4][0] = DR.data_points[3].x; // vert 4
    tv[4][1] = DR.data_points[3].y;
    tv[4][2] = DR.data_points[3].z;
    // midpoints of edges  vert vert  index
    tv[5][0] =  (tv[1][0]+tv[2][0])/2.0; // 12  ix 5 midpoints
    tv[5][1] =  (tv[1][1]+tv[2][1])/2.0;
    tv[5][2] =  (tv[1][2]+tv[2][2])/2.0;
    tv[6][0] = (tv[1][0]+tv[3][0])/2.0; // 13 ix 6
    tv[6][1] = (tv[1][1]+tv[3][1])/2.0;
    tv[6][2] = (tv[1][2]+tv[3][2])/2.0;
    tv[7][0] = (tv[1][0]+tv[4][0])/2.0; // 14 ix 7
    tv[7][1] = (tv[1][1]+tv[4][1])/2.0;
    tv[7][2] = (tv[1][2]+tv[4][2])/2.0;
    tv[8][0] = (tv[2][0]+tv[3][0])/2.0; // 23 ix 8
    tv[8][1] = (tv[2][1]+tv[3][1])/2.0;
    tv[8][2] = (tv[2][2]+tv[3][2])/2.0;
    tv[9][0] = (tv[2][0]+tv[4][0])/2.0; // 24 ix 9
    tv[9][1] = (tv[2][1]+tv[4][1])/2.0;
    tv[9][2] = (tv[2][2]+tv[4][2])/2.0;
    tv[10][0] = (tv[3][0]+tv[4][0])/2.0; // 34 ix 10
    tv[10][1] = (tv[3][1]+tv[4][1])/2.0;
    tv[10][2] = (tv[3][2]+tv[4][2])/2.0;
    tv[11][0] = (tv[1][0]+tv[2][0]+tv[3][0]+tv[4][0])/4.0; // center ix 11
    tv[11][1] = (tv[1][1]+tv[2][1]+tv[3][1]+tv[4][1])/4.0;
    tv[11][2] = (tv[1][2]+tv[2][2]+tv[3][2]+tv[4][2])/4.0;
    ntv = 11;

    System.out.println(" ");
    System.out.println("split into tetrahedrons");
    // tetras  4 triangles each
    tet[1][0] = 1;
    tet[1][1] = 5;
    tet[1][2] = 6;
    tet[1][3] = 7;
    tet[2][0] = 2;
    tet[2][1] = 5;
    tet[2][2] = 8;
    tet[2][3] = 9;
    tet[3][0] = 3;
    tet[3][1] = 6;
    tet[3][2] = 8;
    tet[3][3] = 10;
    tet[4][0] = 4;
    tet[4][1] = 7;
    tet[4][2] = 9;
    tet[4][3] = 10;

    tet[5][0] = 11;
    tet[5][1] = 5;
    tet[5][2] = 6;
    tet[5][3] = 7;
    tet[6][0] = 11;
    tet[6][1] = 5;
    tet[6][2] = 8;
    tet[6][3] = 9;
    tet[7][0] = 11;
    tet[7][1] = 6;
    tet[7][2] = 8;
    tet[7][3] = 10;
    tet[8][0] = 11;
    tet[8][1] = 7;
    tet[8][2] = 9;
    tet[8][3] = 10;

    tet[9][0] = 11;
    tet[9][1] = 5;
    tet[9][2] = 7;
    tet[9][3] = 9;
    tet[10][0] = 11;
    tet[10][1] = 8;
    tet[10][2] = 9;
    tet[10][3] = 10;
    tet[11][0] = 11;
    tet[11][1] = 6;
    tet[11][2] = 7;
    tet[11][3] = 10;
    tet[12][0] = 11;
    tet[12][1] = 5;
    tet[12][2] = 6;
    tet[12][3] = 8;
    
    ntet = 12;
    
    System.out.println(" ");
    System.out.println("vertex x,y,z");
    for(int i=1; i<=ntv; i++)
    {
	System.out.println("tv vert "+i+"  "+tv[i][0]+"  "+tv[i][1]+
			   "  "+tv[i][2]);
    }            

    System.out.println(" ");
    System.out.println("tetrahedron p1 p2 p3 p4");
    for(int i=1; i<=ntet; i++)
    {
	System.out.println("tetrahedron "+i+"  "+tet[i][0]+"  "+tet[i][1]+
			   "  "+tet[i][2]+"  "+tet[i][3]);
    }            
    System.out.println(" ");
    
    // System.out.println("check edge length");
    System.out.println(" ");
    System.out.println("d 1 2 = "+d(1,2));
    System.out.println("d 2 3 = "+d(2,3));
    System.out.println("d 3 4 = "+d(3,4));
    System.out.println("d 4 1 = "+d(4,1));
    System.out.println("d 5 1 = "+d(5,1));
    System.out.println("d 5 2 = "+d(5,2));
    System.out.println("d 6 1 = "+d(6,1));
    System.out.println("d 6 3 = "+d(6,3));
    System.out.println("d 7 1 = "+d(7,1));
    System.out.println("d 7 4 = "+d(7,4));
    System.out.println("d 8 2 = "+d(8,2));
    System.out.println("d 8 3 = "+d(8,3));
    System.out.println("d 9 2 = "+d(9,2));
    System.out.println("d 9 4 = "+d(9,4));
    System.out.println("d 10 3 = "+d(10,3));
    System.out.println("d 10 4 = "+d(10,4));
    System.out.println(" ");
    System.out.println("compute volume");
    volumet = volume(1, 2, 3, 4);
    System.out.println("volume total, check sum of volumes ="+volumet);
		       
    volume[1] = volume(tet[1][0], tet[1][1], tet[1][2], tet[1][3]);
    System.out.println("tet volume[1]="+volume[1]);
    volume[2] = volume(tet[2][0], tet[2][1], tet[2][2], tet[2][3]);
    System.out.println("tet volume[2]="+volume[2]);
    volume[3] = volume(tet[3][0], tet[3][1], tet[3][2], tet[3][3]);
    System.out.println("tet volume[3]="+volume[3]);
    volume[4] = volume(tet[4][0], tet[4][1], tet[4][2], tet[4][3]);
    System.out.println("tet volume[4]="+volume[4]);
    volume[5] = volume(tet[5][0], tet[5][1], tet[5][2], tet[5][3]);
    System.out.println("tet volume[5]="+volume[5]);
    volume[6] = volume(tet[6][0], tet[6][1], tet[6][2], tet[6][3]);
    System.out.println("tet volume[6]="+volume[6]);
    volume[7] = volume(tet[7][0], tet[7][1], tet[7][2], tet[7][3]);
    System.out.println("tet volume[7]="+volume[7]);
    volume[8] = volume(tet[8][0], tet[8][1], tet[8][2], tet[8][3]);
    System.out.println("tet volume[8]="+volume[8]);
    volume[9] = volume(tet[9][0], tet[9][1], tet[9][2], tet[9][3]);
    System.out.println("tet volume[9]="+volume[9]);
    volume[10] = volume(tet[10][0], tet[10][1], tet[10][2], tet[10][3]);
    System.out.println("tet volume[10]="+volume[10]);
    volume[11] = volume(tet[11][0], tet[11][1], tet[11][2], tet[11][3]);
    System.out.println("tet volume[11]="+volume[11]);
    volume[12] = volume(tet[12][0], tet[12][1], tet[12][2], tet[12][3]);
    System.out.println("tet volume[12]="+volume[12]);
    for(int i=1; i<=ntet; i++) { volumetot += volume[i]; }
    System.out.println("volume sum, check sum of volumes ="+volumetot);
    
    System.out.println(" ");

    // set up for datwrite  tet to polys    tv to points
    kp = 0;
    for(int i=1; i<=ntet; i++) // tetra
    {
      for(int k=0; k<4; k++) // 4 faces
      {
        polys[kp] = 3;
	kp++;
	for(int j=k; j<k+3; j++) // triangle for each face
	{
	  jj = j;
	  if(jj>3) jj = jj-4;
	  polys[kp] = tet[i][jj];
	  kp++;
	} // end j
      } // end k
    } // end i 

    // set up for datwrite points
    for(int i=1; i<=ntv; i++)
    {
      for(int j=0; j<3; j++)
      {
	points[i-1][j] = tv[i][j]; // make 0 offset
      }
    }
    // write 4 triangles for each tetrahedron
    new datwrite(fileout, points, ntv, polys, 4*ntet);

    System.out.println("tetra_split1.dat file written");
    System.out.println("tetra_split1.java finished");
  } // end tetra_split1

  double d(int i, int j) // distance between vert i and vert j
  {
    double dist = 0.0;
    dist = Math.sqrt((tv[i][0]-tv[j][0])*(tv[i][0]-tv[j][0])+
		     (tv[i][1]-tv[j][1])*(tv[i][1]-tv[j][1])+
		     (tv[i][2]-tv[j][2])*(tv[i][2]-tv[j][2]));
    return dist;
  }

  double volume(int p1, int p2, int p3, int p4)
  {
    double A[][] = new double[3][3];
    double vol1;
    A[0][0] = tv[p1][0]-tv[p4][0]; 
    A[0][1] = tv[p1][1]-tv[p4][1]; 
    A[0][2] = tv[p1][2]-tv[p4][2]; 
    A[1][0] = tv[p2][0]-tv[p4][0]; 
    A[1][1] = tv[p2][1]-tv[p4][1]; 
    A[1][2] = tv[p2][2]-tv[p4][2]; 
    A[2][0] = tv[p3][0]-tv[p4][0]; 
    A[2][1] = tv[p3][1]-tv[p4][1]; 
    A[2][2] = tv[p3][2]-tv[p4][2];
    vol1 = Math.abs(Matrix.determinant(A)/6.0);
    return vol1;
  } // end volume of tetrahedron
    
  public static void main (String[] args)
  {
      new tetra_split1(args[0],args[1]); // .dat 8 points, .dat out
  }
} // end tetra_split1.java