// tetra_split.java                 check split into smaller tetrahedrons
//                          1
//                          *
// split                   /|\
// 1 2 3 4 center vert    / | \
// 5 6 7 8 outer verts   /  |  \
// 1 2 3 5 outer tet    /   |   \
// 1 2 4 6 outer tet   /   5|    \
// 1 3 4 7 outer tet  /     |     \
// 2 3 4 8 outer te 2*------|------*3
//                  / \     |     / \
// cntr tri        /   \    |    /   \
//                /     \   |8  /     \ 
//               /       \  |  /       \
//              /         \ | /         \
//             /           \|/           \
//           6*------------ *-------------*7
//                          4
//

// needs Matrix.determinant

public class tetra_split
{
  int vert[][] = new int[40][3]; //  triangle faces, 4 per tetra, 1 offset
  int px = 20; // up spare
  int nvert = 4;
  int ncntr = 6;
  double tv[][]   = new double[50][3]; // x,y,z for 4 vertices,
  double points[][] = new double[19][3]; // x,y,z 0 offset
  double mat[][] = new double[nvert][nvert]; // for volume
  double volume = 0.0;
  double fct = 0.54;  // scale outer
  double Pi = Math.PI;
  int polys[] = new int[200]; // at least 4*vert
  int kp = 0;
  int jj = 0;
  int tet[][] = {{1, 2, 3, 4}, // added later 4 vert for each tetra
	 	 {1, 2, 3, 5}, // split
		 {1, 2, 4, 6},
		 {1, 3, 4, 7},
		 {2, 3, 4, 8}};
  double xc, yc, zc;
  double xp, yp, zp;
    
  public tetra_split()
  {
    System.out.println("tetra_split.java running");
    
    vert[1][0] = 1; // tri 1
    vert[1][1] = 2;
    vert[1][2] = 4;
    vert[2][0] = 2; // tri 2
    vert[2][1] = 3;
    vert[2][2] = 4;
    vert[3][0] = 3; // tri 3
    vert[3][1] = 1;
    vert[3][2] = 4;
    vert[4][0] = 1; // tri 4
    vert[4][1] = 2;
    vert[4][2] = 3;

    vert[5][0] = 1; // vert 5
    vert[5][1] = 2;
    vert[5][2] = 5;
    vert[6][0] = 1; 
    vert[6][1] = 3;
    vert[6][2] = 5;
    vert[7][0] = 2;
    vert[7][1] = 3;
    vert[7][2] = 5;

    vert[8][0] = 1; // vert 6
    vert[8][1] = 2;
    vert[8][2] = 6;
    vert[9][0] = 1; 
    vert[9][1] = 4;
    vert[9][2] = 6;
    vert[10][0] = 2;
    vert[10][1] = 4;
    vert[10][2] = 6;

    vert[11][0] = 1; // vert 7
    vert[11][1] = 3;
    vert[11][2] = 7;
    vert[12][0] = 1; 
    vert[12][1] = 4;
    vert[12][2] = 7;
    vert[13][0] = 3;
    vert[13][1] = 4;
    vert[13][2] = 7;

    vert[14][0] = 2; // vert 8
    vert[14][1] = 3;
    vert[14][2] = 8;
    vert[15][0] = 2; 
    vert[15][1] = 4;
    vert[15][2] = 8;
    vert[16][0] = 3;
    vert[16][1] = 4;
    vert[16][2] = 8;

    System.out.println("triangle faces vertex numbers");
    for(int i=1; i<=nvert; i++)
    {
	System.out.println("tri "+i+"  "+vert[i][0]+"  "+vert[i][1]+
			   "  "+vert[i][2]);
    }            
    System.out.println(" ");
    tv[0][0] = 0.0;
    tv[0][1] = 0.0;
    tv[0][2] = 0.0;
    
    tv[1][0] = 0.0; // vert 1 x,y,z up
    tv[1][1] = 0.0;
    tv[1][2] = Math.sqrt(2.0);
    tv[2][0] = 1.0; // vert 2 x,y,z X dir
    tv[2][1] = 0.0;
    tv[2][2] = 0.0;
    tv[3][0] = Math.cos(2.0*Pi/3.0); // vert 3 120 deg
    tv[3][1] = Math.sin(2.0*Pi/3.0);
    tv[3][2] = 0.0;
    tv[4][0] = Math.cos(4.0*Pi/3.0); // vert 4 240 deg
    tv[4][1] = Math.sin(4.0*Pi/3.0);
    tv[4][2] = 0.0;

    System.out.println("vertex x,y,z");
    for(int i=1; i<=nvert; i++)
    {
	System.out.println("vert "+i+"  "+tv[i][0]+"  "+tv[i][1]+
			   "  "+tv[i][2]);
    }            
    System.out.println(" ");
    
    System.out.println("check edge length");
    System.out.println("d 1 2 = "+d(1,2));
    System.out.println("d 1 3 = "+d(1,3));
    System.out.println("d 1 4 = "+d(1,4));
    System.out.println("d 2 3 = "+d(2,3));
    System.out.println("d 2 4 = "+d(2,4));
    System.out.println("d 3 4 = "+d(3,4));
    System.out.println(" ");
    System.out.println("compute area");
    System.out.println("area = "+area(1, 2, 3, 4));
    System.out.println(" ");
    System.out.println("compute volume");
    System.out.println("volume = "+volume(1, 2, 3, 4));
    System.out.println("ck vol = "+(tri_area(2,3,4)*tv[1][2]/3.0));
    System.out.println("ck2 vol = "+volume2(1, 2, 3, 4));
    System.out.println("ck3 vol = "+(d(1,2)*d(1,2)*d(1,2)*Math.sqrt(2.0)/12.0));
    System.out.println(" ");

    System.out.println("split into tetrahedrons");
    System.out.println(" ");

    // 2 3 4 8 
    xc = (tv[2][0]+tv[3][0]+tv[4][0])/3.0;
    yc = (tv[2][1]+tv[3][1]+tv[4][1])/3.0;
    zc = (tv[2][2]+tv[3][2]+tv[4][2])/3.0;
    System.out.println("2 3 4 center "+xc+" "+yc+" "+zc);
    differ(3, 2, px);
    System.out.println("32 diff "+tv[px][0]+" "+tv[px][1]+" "+tv[px][2]);
    differ(4, 2, px+1);
    System.out.println("42 diff "+tv[px+1][0]+" "+tv[px+1][1]+" "+tv[px+1][2]);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("32 42 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("8 vert "+xp+" "+yp+" "+zp);

    differ(2, 3, px);
    differ(2, 4, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("23 24 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("8 vert "+xp+" "+yp+" "+zp);

    differ(2, 3, px);
    differ(3, 4, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("23 34 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("8 vert "+xp+" "+yp+" "+zp);

    differ(4, 2, px);
    differ(4, 3, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("42 43 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("8 vert "+xp+" "+yp+" "+zp);

    differ(2, 4, px);
    differ(4, 3, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("24 43 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("8 vert "+xp+" "+yp+" "+zp);
    tv[8][0] = xp;
    tv[8][1] = yp;
    tv[8][2] = zp;
    System.out.println("check edge length");
    System.out.println("d 2 8 = "+d(2,8));
    System.out.println("d 3 8 = "+d(3,8));
    System.out.println("d 4 8 = "+d(4,8));
    System.out.println(" ");

    // 1 2 3 5 
    xc = (tv[1][0]+tv[2][0]+tv[3][0])/3.0;
    yc = (tv[1][1]+tv[2][1]+tv[3][1])/3.0;
    zc = (tv[1][2]+tv[2][2]+tv[3][2])/3.0;
    System.out.println("1 2 3 center "+xc+" "+yc+" "+zc);
    differ(1, 2, px);
    differ(1, 3, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("12 13 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("5 vert "+xp+" "+yp+" "+zp);

    System.out.println("reverse differ order");
    differ(2, 1, px);
    differ(3, 1, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("5 vert "+xp+" "+yp+" "+zp);

    differ(2, 1, px);
    differ(2, 3, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("21 23 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("5 vert "+xp+" "+yp+" "+zp);

    differ(1, 2, px);
    differ(2, 3, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("12 23 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("5 vert "+xp+" "+yp+" "+zp);
    System.out.println(" ");
    tv[5][0] = xp;
    tv[5][1] = yp;
    tv[5][2] = zp;

    // 1 2 4 6 
    xc = (tv[1][0]+tv[2][0]+tv[4][0])/3.0;
    yc = (tv[1][1]+tv[2][1]+tv[4][1])/3.0;
    zc = (tv[1][2]+tv[2][2]+tv[4][2])/3.0;
    System.out.println("1 2 4 center "+xc+" "+yc+" "+zc);

    differ(1, 2, px);
    differ(2, 4, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("12 24 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc - fct*tv[px+2][0];
    yp = yc - fct*tv[px+2][1];
    zp = zc - fct*tv[px+2][2];
    System.out.println("6 vert "+xp+" "+yp+" "+zp);
    tv[6][0] = xp;
    tv[6][1] = yp;
    tv[6][2] = zp;
    System.out.println("check edge length");
    System.out.println("d 1 6 = "+d(1,6));
    System.out.println("d 2 6 = "+d(2,6));
    System.out.println("d 4 6 = "+d(4,6));
    System.out.println(" ");


    // 1 3 4 7 
    xc = (tv[1][0]+tv[3][0]+tv[4][0])/3.0;
    yc = (tv[1][1]+tv[3][1]+tv[4][1])/3.0;
    zc = (tv[1][2]+tv[3][2]+tv[4][2])/3.0;
    System.out.println("1 3 4 center "+xc+" "+yc+" "+zc);

    differ(1, 3, px);
    differ(3, 4, px+1);
    cross(px, px+1, px+2); // px+2 is temp cross product
    System.out.println("13 34 cross "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    xp = xc + fct*tv[px+2][0];
    yp = yc + fct*tv[px+2][1];
    zp = zc + fct*tv[px+2][2];
    System.out.println("7 vert "+xp+" "+yp+" "+zp);
    tv[7][0] = xp;
    tv[7][1] = yp;
    tv[7][2] = zp;
    System.out.println("check edge length");
    System.out.println("d 1 7 = "+d(1,7));
    System.out.println("d 3 7 = "+d(3,7));
    System.out.println("d 4 7 = "+d(4,7));
    System.out.println(" ");

    System.out.println("check outer length");
    System.out.println("d 5 6 = "+d(5,6));
    System.out.println("d 5 7 = "+d(5,7));
    System.out.println("d 5 8 = "+d(5,8));
    System.out.println("d 6 7 = "+d(6,7));
    System.out.println("d 6 8 = "+d(6,8));
    System.out.println("d 7 8 = "+d(7,8));
    System.out.println(" ");

    System.out.println("check straight edges");
    differ(5, 3, px);
    differ(3, 7, px+1);
    differ(px, px+1, px+2);
    System.out.println("5 3 7 "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    differ(6, 4, px);
    differ(4, 7, px+1);
    differ(px, px+1, px+2);
    System.out.println("6 4 7 "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    differ(5, 2, px);
    differ(2, 8, px+1);
    differ(px, px+1, px+2);
    System.out.println("5 2 8 "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    differ(5, 1, px);
    differ(1, 6, px+1);
    differ(px, px+1, px+2);
    System.out.println("5 1 6 "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    differ(7, 3, px);
    differ(3, 8, px+1);
    differ(px, px+1, px+2);
    System.out.println("7 3 8 "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    differ(8, 3, px);
    System.out.println("diff 83 "+tv[px][0]+" "+tv[px][1]+" "+tv[px][2]);
    differ(3, 7, px+1);
    System.out.println("diff 37 "+tv[px+1][0]+" "+tv[px+1][1]+" "+tv[px+1][2]);
    differ(px, px+1, px+2);
    System.out.println("8 3 7 "+tv[px+2][0]+" "+tv[px+2][1]+" "+tv[px+2][2]);
    System.out.println(" ");
    
    
    System.out.println("split vertex x,y,z");
    for(int i=5; i<=8; i++)
    {
	System.out.println("vert "+i+"  "+tv[i][0]+"  "+tv[i][1]+
			   "  "+tv[i][2]);
    }            
    System.out.println(" ");

    System.out.println("check edge length 1");
    System.out.println("d 1 5 = "+d(1,5));
    System.out.println("d 2 5 = "+d(2,5));
    System.out.println("d 3 5 = "+d(3,5));
    System.out.println("d 1 6 = "+d(1,6));
    System.out.println("d 2 6 = "+d(2,6));
    System.out.println("d 4 6 = "+d(4,6));
    System.out.println(" ");
    System.out.println("d 1 7 = "+d(1,7));
    System.out.println("d 3 7 = "+d(3,7));
    System.out.println("d 4 7 = "+d(4,7));
    System.out.println("d 2 8 = "+d(2,8));
    System.out.println("d 3 8 = "+d(3,8));
    System.out.println("d 4 8 = "+d(4,8));
    System.out.println(" ");
    System.out.println("compute area 1, 2, 3, 5");
    System.out.println("area = "+area(1, 2, 3, 5));
    System.out.println(" ");
    System.out.println("compute volume");
    System.out.println("volume = "+volume(1, 2, 3, 5));
    System.out.println(" ");

    System.out.println("compute area 1, 2, 4, 6");
    System.out.println("area = "+area(1, 2, 4, 6));
    System.out.println(" ");
    System.out.println("compute volume");
    System.out.println("volume = "+volume(1, 2, 4, 6));
    System.out.println("ck vol = "+volck(1, 2, 4, 6));
    System.out.println(" ");

    System.out.println("compute area 1, 3, 4, 7");
    System.out.println("area = "+area(1, 3, 4, 7));
    System.out.println(" ");
    System.out.println("compute volume");
    System.out.println("volume = "+volume(1, 3, 4, 7));
    System.out.println("ck vol = "+volck(1, 3, 4, 7));
    System.out.println(" ");

    kp = 0;
    for(int i=1; i<=16; i++) // all
    {
      polys[kp] = 3;
      kp++;
      for(int j=0; j<3; j++)
      {
	polys[kp] = vert[i][j]; // is 1 offset
	kp++;
      }
    }

    for(int i=1; i<=8; i++)
    {
      for(int j=0; j<3; j++)
      {
	points[i-1][j] = tv[i][j]; // make 0 offset
      }
    }
    new datwrite("tetra_split5.dat", points, 8, polys, 16);

    System.out.println("tetra_split.java finished");
  } // end tetra_split

  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 tri_area(int p1, int p2, int p3)
  {
    double a = d(p1, p2);
    double b = d(p2, p3);
    double c = d(p3, p1);
    double s = (a+b+c)/2.0;
    return Math.sqrt(s*(s-a)*(s-b)*(s-c));
  } // end area of triangle

  double area(int p1, int p2, int p3, int p4)
  {
    double a = 0.0;
    a += tri_area(p1, p2, p3);
    a += tri_area(p2, p3, p4);
    a += tri_area(p3, p4, p1);
    a += tri_area(p4, p1, p2);
    return a;
  } // end area of tetrahedron

  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

  double volume2(int p1, int p2, int p3, int p4)
  {
    double A[][] = new double[4][4];
    double vol1;
    A[0][0] = tv[p1][0]; 
    A[0][1] = tv[p1][1]; 
    A[0][2] = tv[p1][2]; 
    A[0][3] = 1.0; 
    A[1][0] = tv[p2][0]; 
    A[1][1] = tv[p2][1]; 
    A[1][2] = tv[p2][2]; 
    A[1][3] = 1.0; 
    A[2][0] = tv[p3][0]; 
    A[2][1] = tv[p3][1]; 
    A[2][2] = tv[p3][2];
    A[2][3] = 1.0; 
    A[3][0] = tv[p4][0]; 
    A[3][1] = tv[p4][1]; 
    A[3][2] = tv[p4][2];
    A[3][3] = 1.0; 
    vol1 = Math.abs(Matrix.determinant(A)/6.0);
    return vol1;
  } // end volume of tetrahedron


  void cross(int p1, int p2, int p3) // p3 new tv
  {
    tv[p3][0] = tv[p1][1]*tv[p2][2] - tv[p1][2]*tv[p2][1];
    tv[p3][1] = tv[p1][2]*tv[p2][0] - tv[p1][0]*tv[p2][2];
    tv[p3][2] = tv[p1][0]*tv[p2][1] - tv[p1][1]*tv[p2][0];
		// a.dy*b.dz - a.dz*b.dy,
		// a.dz*b.dx - a.dx*b.dz,
		// a.dx*b.dy - a.dy*b.dx);
  } // end  cross

  void differ(int p1, int p2, int p3)
  {
    tv[p3][0] = tv[p1][0]-tv[p2][0];
    tv[p3][1] = tv[p1][1]-tv[p2][1];
    tv[p3][2] = tv[p1][2]-tv[p2][2];
  } // end differ
    
  double dot(int p1, int p2)
  {
    return tv[p1][0]*tv[p2][0] +
           tv[p1][1]*tv[p2][1] +
	   tv[p1][2]*tv[p2][2];
           // a.dx*b.dx + a.dy*b.dy + a.dz*b.dz;
  } // end dot

  double volck(int p1, int p2, int p3, int p4)
  {
    // dot((p1-p4),(cross((p2-p4),(p3-p4))));
    //      px             px+1    px+2
    //                         px+3
    differ(p1, p4, px);
    differ(p2, p4, px+1);
    differ(p3, p4, px+2);
    cross(px+1, px+2, px+3);
    return dot(px, px+3)/6.0;
  } // end volck
    
  public static void main (String[] args)
  {
    new tetra_split();
  }
} // end tetra_split.java