// test_tet_split.java    compile    javac -cp . tetra_split.java 
//                        compile    javac -cp . datread.java
//                        compile    javac -cp . datwrite.java
//                        compile    javac -cp . test_tet_split.java
//                        execute    java  -cp . test_tet_split
// needs datread.class datwrite.class and input file  xxx.tet

public class test_tet_split
{         
  int ntet[] = new int[1];                // number of tetra passed and updated
  int tpts[][] = new int[2000][4];        // four point numbers 0..ntet[0]-1
  int ncoord[] = new int[1];              // number of coord passed and updated
  double tcoord[][] = new double[500][3]; // three coordinates x,y,z 
  int stet;

  public test_tet_split(String tet_file, String out_file)
  {
    int pts, pt, k;
    datread DR; // = new datread(tet_file); // defines data structures

    System.out.println("test_tet_split.java running");
    System.out.println("splitting file "+tet_file);
    System.out.println("writing file "+out_file);

    // read vertices into tcoord[vert][0] [1] [2] x, y, z
    System.out.println("datread.java reading "+tet_file);
    DR = new datread(tet_file); // defines data structures

    System.out.println("num_points="+DR.num_points+
                       ", num_polys="+DR.num_polys);
    System.out.println("print vertices");
    for(int i=0; i<DR.num_points; i++) // vertices
    {
      System.out.println(+(i+1)+"  "+DR.data_points[i].x+
                         "  "+DR.data_points[i].y+"  "+DR.data_points[i].z);
      tcoord[i][0] = DR.data_points[i].x;
      tcoord[i][1] = DR.data_points[i].y;
      tcoord[i][2] = DR.data_points[i].z;

    } // end print vertices

    // read polygons into tpts[poly][0] [1] [2] [3] -1 for index tcoord[][]
    System.out.println("print polygons");
    k = 0; // polygons
    for(int i=0; i<DR.num_polys; i++)
    {
      pts = DR.data_polys[k++];
      System.out.print(+pts); // no cr
      for(int j=0; j<pts; j++)
      {
        pt = DR.data_polys[k++];
        tpts[i][j] = pt-1;
        System.out.print(" "+ pt);
      }
      System.out.println(" ");
    } // end print polygons


    ntet[0] = DR.num_polys;
    ncoord[0] = DR.num_points; // will grow
    stet = 0; // index of tetra to split
    for(int stet=0; stet<DR.num_polys; stet++)
      new tetra_split(stet, ntet, tpts,  // four point numbers
                      ncoord, tcoord);   // three coordinates x,y,z 

    System.out.println(" ");
    System.out.println("results after tetra_split ");
    System.out.println(" ");
    System.out.println("tet, p1, p2, p3, p4   index");
    for(int i=0; i<ntet[0]; i++)
	System.out.println(" "+i+",  "+tpts[i][0]+", "+tpts[i][1]+
			   ", "+tpts[i][2]+", "+tpts[i][3]);

    System.out.println(" ");
    System.out.println("pt, x, y, z");
    for(int i=0; i<ncoord[0]; i++)
	System.out.println(" "+i+", "+tcoord[i][0]+", "+tcoord[i][1]+
			   ", "+tcoord[i][2]);
    System.out.println(" ");
    write_dat(out_file, ntet, tpts,        // four point numbers
              ncoord, tcoord);   // three coordinates x,y,z 

    System.out.println(" ");
    System.out.println("check volume ");
    check_volume(ntet, tpts,        // four point numbers
                 ncoord, tcoord);   // three coordinates x,y,z 

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


  void write_dat(String filename, int ntet[], int tpts[][], // four pt numbers
	         int ncoord[], double tcoord[][]) // three coordinates x,y,z
  {
    int num_points = 0;
    int num_polys = 0;;
    double data_points[][] = new double[500][3];
    int data_verts[] = new int[2000];
    int k = 0; // data_verts has entries  4  p1 p2 p3 p4 points +1 of coord

    for(int i=0; i<ncoord[0]; i++)
      for(int j=0; j<3; j++)
	data_points[i][j] = tcoord[i][j];
    num_points = ncoord[0];

    for(int i=0; i<ntet[0]; i++)
    {
      data_verts[k++] = 4;
      data_verts[k++] = tpts[i][0]+1; 
      data_verts[k++] = tpts[i][1]+1; 
      data_verts[k++] = tpts[i][2]+1; 
      data_verts[k++] = tpts[i][3]+1; 
    }
    num_polys = ntet[0]; // 4 4 pts  for each tetrahedron
    new datwrite(filename, data_points, num_points,
                           data_verts, num_polys);
  } // end write_dat

  void check_volume(int ntet[], int tpts[][],      // four point numbers
	          int ncoord[], double tcoord[][]) // three coordinates x,y,z
  {
    double B[][] = new double[4][4];
    double vol1234, vol_sum, vol;
    int nt;

    for(int i=0; i<4; i++) B[i][3] = 1.0;

    vol_sum = 0.0;
    for(int k=0; k<ntet[0]; k++)
    {
      System.out.print("k="+k+",");
      for(int i=0; i<4; i++) B[i][3] = 1.0;
      for(int i=0; i<4; i++)
      {
	nt = tpts[k][i];
        System.out.print("  "+nt);
        for(int j=0; j<3; j++) B[i][j] = tcoord[nt][j];
      } // end i
      System.out.println(" ");

      vol = Math.abs(determinant(B)/6.0);
      System.out.println("k="+k+", vol="+vol);
      vol_sum = vol_sum + vol;
    } // end k
    System.out.println("vol_sum="+vol_sum);

  } // end check_volume

  double determinant(double A[][])
  {
    int n=A.length;
    double D = 1.0;                 // determinant
    double B[][]=new double[n][n];  // working matrix
    int row[]=new int[n];           // row interchange indicies
    int hold , I_pivot;             // pivot indicies
    double pivot;                   // pivot element value
    double abs_pivot;

    if(A[0].length!=n)
    {
      System.out.println("Error in Matrix.determinant, inconsistent array sizes.");
    }
    // build working matrix
    for(int i=0; i<n; i++)
      for(int j=0; j<n; j++)
        B[i][j]=A[i][j];
    // set up row interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k]= k;
    }
    // begin main reduction loop
    for(int k=0; k<n-1; k++)
    {
      // find largest element for pivot
      pivot = B[row[k]][k];
      abs_pivot = Math.abs(pivot);
      I_pivot = k;
      for(int i=k; i<n; i++)
      {
        if( Math.abs(B[row[i]][k]) > abs_pivot )
        {
          I_pivot = i;
          pivot = B[row[i]][k];
          abs_pivot = Math.abs(pivot);
        }
      }
      // have pivot, interchange row indicies
      if(I_pivot != k)
      {
        hold = row[k];
        row[k] = row[I_pivot];
        row[I_pivot] = hold;
        D = - D;
      }
      // check for near singular
      if(abs_pivot < 1.0E-20)
      {
        return abs_pivot;
      }
      else
      {
        D = D * pivot;
        // reduce about pivot
        for(int j=k+1; j<n; j++)
        {
          B[row[k]][j] = B[row[k]][j] / B[row[k]][k];
        }
        //  inner reduction loop
        for(int i=0; i<n; i++)
        {
          if(i != k)
          {
            for(int j=k+1; j<n; j++)
            {
              B[row[i]][j] = B[row[i]][j] - B[row[i]][k]* B[row[k]][j];
            }
          }
        }
      }
      //  finished inner reduction
    }
    // end of main reduction loop
    return D * B[row[n-1]][n-1];
  } // end determinant

  public static void main (String[] args)
  {
    // must have input  .tet  file
    // creates output   .tet  file
    new test_tet_split(args[0], args[1]);
  }
} // end test_tet_split.java