// test_eqn33.java   compile    javac -cp . test_eqn33.java
//                   execute    java  -cp . test_eqn33
// solve 2*U(x,y,z)+3*V(x,y,z)+4*W(x,y,z) = C(x,y,z)
// solve 4*U(x,y,z)+2*V(x,y,z)+3*W(x,y,z) = D(x,y,z)
// solve 3*U(x,y,z)+4*V(x,y,z)+2*W(x,y,z) = E(x,y,z)
// for values of U(x,y,z), V(x,y,z), W(x,y,z) at points in a grid
// U(x,y,z) = 1+x*x*x-y+z   V(x,y,z) = x*x+y*y+z*z  W(x,y,z) = 2-x+2y+z*z
// unknown functions
// C, D, E values given
// needs set of 3 equations for 3 variables at each grid
// needs simeq.java compiled


public class test_eqn33
{
  int nx = 2; // grid
  int ny = 3;
  int nz = 2;
  int m, mm, ii, jj, kk;
  int row = 0; // typically S(i,j,k)
  int grid = nx*ny*nz; // points in grid
  int nvar = 3; // U(x,y,z)  V(x,y,z)  W(x,y,z)
  int neqn = nvar*grid; // nvar unknowns at each grid point
  double xmin = 1.0;
  double xmax = 3.0;
  double ymin = 2.0;
  double ymax = 4.0;
  double zmin = 3.0;
  double zmax = 5.0;
  double x, y, z, dx, dy, dz, val;
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double A[][] = new double[neqn][neqn]; // A*X=Y  
  double X[] = new double[neqn];  // unknown  
  double Y[] = new double[neqn];  // known from C(x,y,z),D(x,y,z),E(x,y,z)  
  double err, max_err;
    
  public test_eqn33()
  {
    System.out.println("test_eqn33.java running");
    System.out.println("grid nx="+nx+"  ny="+ny+"   nz="+nz);
    System.out.println("grid="+grid+"  nvar="+nvar+"   neqn="+neqn);
    System.out.println("xmin="+xmin+"  xmax="+xmax);
    System.out.println("ymin="+ymin+"  ymax="+ymax);
    System.out.println("zmin="+zmin+"  zmax="+zmax);
    System.out.println("x grid xg[]");
    xg[0] = xmin;
    dx = (xmax-xmin)/(double)(nx-1);
    System.out.println("dx="+dx+"  xg[0]="+xmin);
    for(int i=1; i<nx; i++)
    {
      xg[i] = xg[i-1]+dx; 
      System.out.println("  xg["+i+"]="+xg[i]);
    } // end i
    
    System.out.println(" ");
    System.out.println("y grid yg[]");
    yg[0] = ymin;
    dy = (ymax-ymin)/(double)(ny-1);
    System.out.println("dy="+dy+"  yg[0]="+yg[0]);
    for(int j=1; j<ny; j++)
    {
      yg[j] = yg[j-1]+dy; 
      System.out.println("  yg["+j+"]="+yg[j]);
    } // end i
    
    System.out.println(" ");
    System.out.println("z grid zg[]");
    zg[0] = zmin;
    dz = (zmax-zmin)/(double)(nz-1);
    System.out.println("dz="+dz+"  zg[0]="+zg[0]);
    for(int k=1; k<nz; k++)
    {
      zg[k] = zg[k-1]+dz; 
      System.out.println("  zg["+k+"]="+zg[k]);
    } // end k
    System.out.println(" ");
    System.out.println("build A and Y ");
    for(int i=0; i<neqn; i++)
    {
      for(int j=0; j<neqn; j++)
      {
	A[i][j] = 0.0;
      } // end j
    } // end i
    m = 0;
    for(int i=0; i<nx; i++)
    {
      x = xg[i]; 
      for(int j=0; j<ny; j++)
      {
        y = yg[j];
        for(int k=0; k<nz; k++)
        {
          z = zg[k];
	  row = S(i,j,k);
	  Y[row] = C(x,y,z);
	  Y[row+grid] = D(x,y,z);
	  Y[row+2*grid] = E(x,y,z);
          System.out.println("U("+x+","+y+","+z+")="+U(x,y,z)+"  V(x,y,z)="+
			     V(x,y,z)+"  W(x,y,z)="+W(x,y,z));

          System.out.println(" at i="+i+",j="+j+",k="+k+" Y["+row+"]="+Y[row]);
          System.out.println("Y["+(row+grid)+"]="+Y[row+grid]);
          System.out.println("Y["+(row+2*grid)+"]="+Y[row+2*grid]);
 	  A[row][m] = 2.0;   // 2*U(x,y,z)
	  A[row][m+1] = 3.0; // 3*V(x,y,z)
	  A[row][m+2] = 4.0; // 4*W(x,y,z)
 	  A[row+grid][m] = 4.0;   // 3*U(x,y,z)  sum is C(x,y,z)
          A[row+grid][m+1] = 2.0; // 2*V(x,y,z)  sum is D(x,y,z)
          A[row+grid][m+2] = 3.0; // 2*V(x,y,z)  sum is E(x,y,z)
	  A[row+2*grid][m] = 3.0;   // 3*U(x,y,z)  sum is C(x,y,z)
          A[row+2*grid][m+1] = 4.0; // 2*V(x,y,z)  sum is D(x,y,z)
          A[row+2*grid][m+2] = 2.0; // 2*V(x,y,z)  sum is E(x,y,z)
	  m = m+nvar;
	} // end k
      } // end j
    } // end i
    System.out.println(" ");
    printA();
    System.out.println("run simeq ");
    new simeq(A,Y,X);
    System.out.println("X results from simeq");
    max_err = 0.0;
    m = 0; // m==0 first variable  m==2 second variable
    for(int row=0; row<neqn; row++)
    {
      mm = row/nvar;
      ii = mm/(ny*nz);
      jj = (mm-ii*ny*nz)/nz;
      kk = (mm - ii*ny*nz - jj*nz);
      System.out.println("ii="+ii+"  jj="+jj+"   kk="+kk+"   m="+m);
      if(m==0) { val = U(xg[ii],yg[jj],zg[kk]); }
      if(m==1) { val = V(xg[ii],yg[jj],zg[kk]); }
      if(m==2) { val = W(xg[ii],yg[jj],zg[kk]); }
      err = Math.abs(X[row]-val);
      if(err>max_err) { max_err=err; }
      System.out.println("X["+row+"]="+X[row]+"  err="+err);
      m = m + 1;
      if(m==nvar) { m=0; }
    } // end row    
    System.out.println(" ");
    System.out.println("max_err="+max_err);
    System.out.println("test_eqn33.java finished");
  } // end test_eqn33

  double U(double x, double y, double z) // for testing and generating matrix
  {
    return 1.0+x*x*x-y+z;
  } // end U

  double V(double x, double y, double z)
  {
    return x*x+y*y+z*z;
  } // end V

  double W(double x, double y, double z)
  {
    return 2.0-x+2.0*y+z*z;
  } // end W

  double C(double x, double y, double z)
  {
    return 2.0*U(x,y,z)+3.0*V(x,y,z)+4.0*W(x,y,z);
  } // end C

  double D(double x, double y, double z)
  {
    return 4.0*U(x,y,z)+2.0*V(x,y,z)+3.0*W(x,y,z);
  } // end D

  double E(double x, double y, double z)
  {
    return 3.0*U(x,y,z)+4.0*V(x,y,z)+2.0*W(x,y,z);
  } // end E
      
  int S(int i, int j, int k) // index to one grid
  {
    return nz*ny*i + nz*j + k;
  } // end S

  void  printA() // non zero values of A
  {
    System.out.println("non zero A");
    for(int i=0; i<neqn; i++)
    {
      for(int j=0; j<neqn; j++)
      {
	if(Math.abs(A[i][j]) > 0.0001)
	{
          System.out.println("A["+i+"]["+j+"]="+A[i][j]);
	} // end if
      } // end j
    } // end i
  } // end printA

  public static void main (String[] args)
  {
    new test_eqn33();
  }
} // end test_eqn33.java