// pde23a_eq.java  3D uniform grid boundary value PDE
//                known solution to test method
//     U(x,y,z) = x*x*x+y*y*y+z*z*z
// 
//     differential equation to solve
//     Uxx + Uyy + Uyzz = c(x,y,z)
//     Ux(x,y,z) = 3 x*x
//     Uy(x,y,z) = 3 y*y
//     Uz(x,y,z) = 3 z*z
//     Uxx(x,y,z) = 6 x
//     Uyy(x,y,z) = 6 y
//     Uzz(x,y,z) = 6 z
//     c(x,y,z) = 6 x + 6 y + 6 z
// 
// uniform grid on cube -1 to 1, -1 to 1, -1 to 1                    
// set up and solve system of (nx-2)*(ny-2)*(nz-2) linear equations
//
// uses nuderiv.java
// uses simeq.java
import java.text.*;
import java.io.*;

class pde23a_eq
{
  int nx = 5;          // includes boundary 
  int ny = 5;          // includes boundary 
  int nz = 5;          // includes boundary 
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double zg[] = new double[nz];
  double cxx[] = new double[nx];
  double cyy[] = new double[ny];
  double czz[] = new double[nz];
  // nx-2 by ny-2 by n-2  internal, non boundary cells 
  int neqn = (nx-2)*(ny-2)*(nz-2);
  int cs = neqn;
  double us[] = new double[neqn];               // solution being computed 
  double ut[][] = new double[neqn][neqn+1];     // equations incl RHS 
  double xmin = -1.0;
  double xmax = 1.0;
  double ymin = -1.0;
  double ymax = 1.0;
  double zmin = -1.0;
  double zmax = 1.0;
  double hx, hy, hz;
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     
  double time_start;
  double time_now;
  double time_last;
  DecimalFormat fe = new DecimalFormat("0.00E00");
  DecimalFormat f6 = new DecimalFormat("0.00000");

  pde23a_eq()
  {
    System.out.println("pde23a_eq.java running");
    System.out.println("differential equation to solve");
    System.out.println("d^2U/dx^2 + d^2U/dy^2 + d^2U/dz^2  = c(x,y,z)");
    System.out.println("Ux(x,y,z) = 3 x*x   for testing");
    System.out.println("Uy(x,y,z) = 3 y*y");
    System.out.println("Uz(x,y,z) = 3 z*z");
    System.out.println("Uxx(x,y,z) = 6 x");
    System.out.println("Uyy(x,y,z) = 6 y");
    System.out.println("Uzz(x,y,z) = 6 z");
    System.out.println("c(x,y,z) = 6 x + 6 y + 6 z  given ");
    System.out.println("uniform grid on cube -1,1 to -1,1 to -1,1");
    System.out.println("known Solution, for testing method");
    System.out.println("U(x,y,z) = x*x*x + y*y*y + z*z*z");
    System.out.println(" and U(x,y,z) in program is for boundary values ");
    System.out.println(" ");

    time_start = (double)System.currentTimeMillis()/1000.0;
    time_last = time_start;
    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);
    hz = (zmax-zmin)/(double)(nz-1);

    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin+(double)i*hx;
      System.out.println("xg("+i+")="+xg[i]);
    } // end i
    for(int j=0; j<ny; j++)
    {
      yg[j] = ymin+(double)j*hy;
      System.out.println("yg("+j+")="+yg[j]);
    } // end j
    for(int k=0; k<nz; k++)
    {
      zg[k] = zmin+(double)k*hz;
      System.out.println("zg("+k+")="+zg[k]);
    } // end k

    System.out.println("xmin="+f6.format(xg[0])+
		       ", xmax="+f6.format(xg[nx-1])+
                       ", nx="+nx+", hx="+f6.format(hx)); 
    System.out.println("ymin="+f6.format(yg[0])+
		       ", ymax="+f6.format(yg[ny-1])+
                       ", ny="+ny+", hy="+f6.format(hy)); 
    System.out.println("zmin="+f6.format(zg[0])+
		       ", zmax="+f6.format(zg[nz-1])+
                       ", nz="+nz+", hz="+f6.format(hz)); 
    System.out.println("neqn="+neqn);
    System.out.println("matrix ut[neqn][neqn+1] "+(neqn*(neqn+1)));
    System.out.println("u(0.5,0.5,0.5)="+
                       f6.format(u(0.5,0.5,0.5)));
    System.out.println("c(0,5,0.5,0.5)="+
                       f6.format(c(0.5,0.5,0.5)));
    System.out.println("S(1,1,1)="+S(1,1,1));
    System.out.println("S(1,2,3)="+S(1,2,3));
    System.out.println("S(nx-1,ny-1,nz-1)="+S(nx-1,ny-1,nz-1));

    init_matrix();
    System.out.println("matrix initialized");
    System.out.println("initial matrix, left side, upper ");
    // printB();
    check_rows(); 
    new simeq(neqn, ut, us);
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time to compute solution = "+(time_now-time_last)+
                       " seconds ");
    time_last = time_now;
    check_soln();
    print();

    check();
    avgerror = error/(double)neqn;
    System.out.println("total error="+fe.format(error)+
                       ", average error="+fe.format(avgerror)+
                       ", max error="+fe.format(maxerror));
    write_soln3(us);
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("total time = "+(time_now-time_start)+
                       " seconds ");
    System.out.println("finished pde23a_eq.java");
  } // end pde23a_eq


  double u(double x, double y, double z) // for boundary and optional check 
  {
    return x*x*x+y*y*y+z*z*z;
  } // end u

  double c(double x, double y, double z)
  // RHS of differential equation to solve 
  {
    return 6.0*x + 6.0*y + 6.0*z;
  } // end c

  int S(int i, int j, int k) // zero based index includes bound
  {                           // this is row of ut[]
    return (i-1)*(ny-2)*(nz-2)+(j-1)*(nz-2)+(k-1);
    // index into us and ut, RHS, no boundary
  } // end S

  void check() // typically not known, instrumentation 
  {
    double uexact, err;
  
    System.out.println("check");
    error = 0.0;
    maxerror = 0.0;
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      { 
        for(int k=1; k<nz-1; k++)
        { 
	  uexact = u(xg[i], yg[j], zg[k]);
          err = Math.abs(us[S(i,j,k)]-uexact);
          error = error + err;
          if(err>maxerror) maxerror=err;
	} // end k
      } // end j
    } // end i
    System.out.println("check against known solution, maxerror="+maxerror);
    System.out.println(" ");
  } // end check 

  void check_rows() // only for(checking 
  {	// row is in ut solution coordinates 
    int row, col;
    double sum;

    System.out.println("check_rows");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  row = S(i,j,k);
          sum = 0.0;
          for(int ii=1; ii<nx-1; ii++)
          {
            for(int jj=1; jj<ny-1; jj++)
	    {
              for(int kk=1; kk<nz-1; kk++)
	      {
		col = S(ii,jj,kk);
                sum = sum + ut[row][col]*u(xg[ii],yg[jj],zg[kk]);
		// u is check solution
	      } // end kk
	    } // jj 
          } // ii 
          sum = sum - ut[row][cs];
          if(Math.abs(sum)>0.001)
          {
            System.out.println("row i="+i+", j="+j+", is bad err="+sum);
          }
	} // end k
      } // end j 
    } // end i 
  } // end check_rows 

  void print() // typically not known 
  {
    double error = 0.0;
    double max_error = 0.0;
    double avg_error = 0.0;
  
    System.out.println("exact solution u, computed us, error");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  error = u(xg[i],yg[j],zg[k]) - us[S(i,j,k)];
          avg_error = avg_error + Math.abs(error);
          if(Math.abs(error)>max_error) max_error = Math.abs(error);

	  System.out.println("xg["+i+"]="+f6.format(xg[i])+
	      	             ", yg["+j+"]="+f6.format(yg[j])+
	      	             ", zg["+k+"]="+f6.format(zg[k])+
			     ", u="+f6.format(u(xg[i],yg[j],zg[k]))+
                             ", us="+f6.format(us[S(i,j,k)])+
                             ", err="+fe.format(error));
	} // end k
      } // end j
    } // end i
    System.out.println("max_error="+max_error+", avg_error="+
                       avg_error/(double)neqn);
    System.out.println(" ");
  } // end print 

  void printB()
  {
    System.out.println("matrix B includes RHS");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
          for(int ii=1; ii<nx-1; ii++)
          {
            for(int jj=1; jj<ny-1; jj++)
	    {
              for(int kk=1; kk<nz-1; kk++)
	      {
                System.out.println("i="+i+", j="+j+", k="+k+
				   " ii="+ii+", jj="+jj+", kk="+kk+
				   ", ut="+ut[S(i,j,k)][S(ii,jj,jj)]);
	      } // end kk
	    } // end jj 
          } // end ii 
          System.out.println("RHS=                         "+ut[S(i,j,k)][cs]);
	} // end k
      } // end j 
    } // end i 
    System.out.println("");
    System.out.println("");
  } //  end printB 

  void init_matrix()
  {
    int row, iii, jjj, kkk, ijj, ikk;

    // cs is RHS, constant column subscript 
    // zero internal cells 
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
          for(int ii=1; ii<nx-1; ii++)
          {
            for(int jj=1; jj<ny-1; jj++)
	    {
              for(int kk=1; kk<nz-1; kk++)
              {
		  ut[S(i,j,k)][S(ii,jj,kk)] = 0.0;
	      } // end kk
            } // end jj
          } // end ii
          ut[S(i,j,k)][cs] = 0.0;
	} // end k
      } // end j
    } // end i

    System.out.println("internal cells zeroed ");
    for(int i=1; i<nx-1; i++) // inside boundary 
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int k=1; k<nz-1; k++)
        {
	  row = S(i,j,k); // equation 
          // for(each term 

          // d^2 U(x,y,z)/dx^2 
          new nuderiv(2,nx,i,xg,cxx);
          for(int ii=0; ii<nx; ii++)
          {
            if(ii==0 || ii==nx-1)
	    {
	      ut[row][cs] = ut[row][cs] - cxx[ii]*u(xg[ii],yg[j],zg[k]);
            }
            else
	    {
	      iii = S(ii,j,k);
              ut[row][iii] = ut[row][iii] + cxx[ii];
	    }
          } // end ii 

          //  d^2 U(x,y,z)/dy^2 
          new nuderiv(2,ny,j,yg,cyy);
          for(int jj=0; jj<ny; jj++)
          {
            if(jj==0 || jj==ny-1)
	    {
	      ut[row][cs] = ut[row][cs] - cyy[jj]*u(xg[i],yg[jj],zg[k]);
            }
            else
	    {
	      ijj = S(i,jj,k);
              ut[row][ijj] = ut[row][ijj] + cyy[jj];
	    }
          } // jj

          // d^2 U(x,y,z)/dz^2 
          new nuderiv(2,nz,k,zg,czz);
          for(int kk=0; kk<nz; kk++)
          {
            if(kk==0 || kk==nz-1)
	    {
	      ut[row][cs] = ut[row][cs] - czz[kk]*u(xg[i],yg[j],zg[kk]);
            }
            else
	    {
	      ikk = S(i,j,kk);
              ut[row][ikk] = ut[row][ikk] + czz[kk];
	    }
          } // end kk

          // RHS constant 
	  ut[row][cs] = ut[row][cs] + c(xg[i],yg[j],zg[k]);
          // System.out.println("i="+i+" j="+j+" k="+k+" row="+row);
	  // System.out.println("ut[row][cs]="+ut[row][cs]);
          // end terms
	} // end k
      } // j 
    } // i 
  } // end init_matrix 

  void check_soln()  // check when solution unknown 
  {
    double smaxerr, err, x, y, z, ux, uy, uz, uxx, uxy, uxz, uyy, uyz, uzz;
    // through dof equations 
    // Uxx + Uyy + Uzz = c(x,y,z) 
    System.out.println("check_soln against PDE");
    smaxerr = 0.0;

    for(int i=1; i<nx-1; i++)  
    {
      new nuderiv(2,nx,i,xg,cxx);
      x = xg[i];
      for(int j=1; j<ny-1; j++)
      {
        new nuderiv(2,ny,j,yg,cyy);
        y = yg[j];
        for(int k=1; k<nz-1; k++)
        {
          new nuderiv(2,nz,k,zg,czz);
          z = zg[k];
	
	  // compute numerical derivative at (i,j,k) 
	  uxx = 0.0;
          for(int ii=0; ii<nx; ii++)
          {
            if(ii==0 || ii==nx-1)
	    {
	      uxx = uxx + cxx[k]*u(xg[ii],yg[j],zg[k]); // boundary 
	    }
            else
	    {
              uxx = uxx + cxx[k]*us[S(ii,j,k)];
            } // end if
          } // end ii
          uyy = 0.0;
          for(int jj=0; jj<ny; jj++)
          {
            if(jj==0 || jj==ny-1)
	    {
              uyy = uyy + cyy[jj]*u(xg[i],yg[jj],zg[k]); // boundary 
	    }
            else
	    {
		uyy = uyy + cyy[k]*us[S(i,jj,k)];
	    } // end if
          } // end jj
          uzz = 0.0;
          for(int kk=0; kk<nz; kk++)
          {
            if(kk==0 || kk==nz-1)
	    {
              uzz = uzz + czz[kk]*u(xg[i],yg[i],zg[kk]); // boundary 
	    }
            else
	    {
		uzz = uzz + czz[k]*us[S(i,j,kk)];
	    } // end if
          } // end kk

	  // coefficients included above
	  err = -c(x,y,z); // add other side of equation terms
	  err = err + uxx + uyy + uzz;
          if(Math.abs(err) > smaxerr) smaxerr = Math.abs(err);
	} // k z
      } // j y 
    } // i x 

    System.out.println("check soln against PDE max error="+smaxerr);
  } // end Check_Soln 

  // write_soln3  for gnuplot 
  void  write_soln3(double Ux[])
  {
    try
    {
      FileWriter fout = new FileWriter("pde23a_eq_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde23a_eq_java.dat");

      for(int i=0; i<nx; i++)
      {
        for(int j=0; j<ny; j++)
        {
          for(int k=0; k<nz; k++)
          {
            if(i==0 || i==nx-1 || j==0 || j==ny-1 || k==0 || k==nz-1)
            {
              fileout.write(xg[i]+" "+yg[j]+" "+zg[k]+" "+
			    u(xg[i],yg[j],zg[k])+"\n"); // boundary 
            }
            else
            { 
              fileout.write(xg[i]+" "+yg[j]+" "+zg[k]+" "+
			    Ux[S(i,j,k)]+"\n"); // computed solution 
            } // end if
	  } // end k
          fileout.write("\n"); // blank line 
        } // end j
        fileout.write("\n"); // blank line 
      } // end i 
      fileout.close();
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde23a_eq_java.dat  not opened");
    }
    catch(Exception exception)
    {
      System.out.println("pde23a_eq_java.dat  not opened");
    }
    System.out.println("finished writing pde23a_eq_java.dat");
  }  // end write_soln3

  public static void main (String[] args)
  {
    new pde23a_eq();
  }
} // end class pde23a_eq