// pdenu34_eq.java  3D non uniform grid boundary value PDE  non iterative 
//                 4th order, exponential, log step size test
//                 known solution to test method                         
//                 u(x,y,z) = exp(x*y*z)
//
//                 differential equation to solve                       
// uxxxx(x,y,z)+2*uyyyy(x,y,z)+3*uzzzz(x,y,z)+4*uxyz(x,y,z) = f(x,y,z)
// f(x,y,z) = (y*y*y*y*z*z*z*z+
//             2*x*x*x*x*z*z*z*z+
//             3*x*x*x*x*y*y*y*y+
//             4+12*z*x*y+
//             4*y*y*z*z*x*x)*exp(x*y*z)
//
// non uniform grid
// Try log spacing nuderiv  0.0, 0.25, 0.5, 0.7, 0.85, 0.95, 1.0
// step = 0.1 then 0.01
// xg[0] = 2.0
// until xg[i-1]-step<0.0  xg[i++] = xg[i-1]-step;  step=2*step  
// xg[i++] = 0.0
//                                                                   
// set up and solve system of linear equations                       
//                                                                   
// create two dimensional array with                                 
// (nx-2)*(ny-2)*(nz-2) rows, one for each equation                         

// uses nuderiv.java
// uses simeq.java
import java.text.*;

class pdenu34_eq
{
  double xg[] = {0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0};
  double yg[] = {0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0};
  double zg[] = {0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1.0};
  int nx = xg.length;          // includes boundary 
  int ny = yg.length;          // includes boundary 
  int nz = zg.length;          // includes boundary 
                      // nx-2 by ny-2 by nz-2  internal, non boundary cells 
  double cx[] = new double[nx];    // nuderiv coefficients 
  double cy[] = new double[ny];
  double cz[] = new double[nz];
  double cxyz[][][] = new double[nx][ny][nz];
  double cxxxx[] = new double[nx];
  double ccxxxx[] = new double[nx]; // for checking
  double cyyyy[] = new double[ny];
  double czzzz[] = new double[nz];
  double coefdxxxx, coefdyyyy, coefdzzzz, coefdxdydz;
  int nxyz = (nx-2)*(ny-2)*(nz-2);
  double us[] = new double[nxyz];              // solution being computed  
  double ut[][] = new double[nxyz][nxyz+1];    // temporary equations 
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     
  DecimalFormat fe = new DecimalFormat("0.00E00");
  DecimalFormat f6 = new DecimalFormat("0.00000");

  pdenu34_eq()
  {

    System.out.println("pdenu34_eq.java running");
    System.out.println("differential equation to solve");
    System.out.println("uxxxx(x,y,z)+2*uyyyy(x,y,z)+3*uzzzz(x,y,z)+4*uxyz(x,y,z) = f(x,y,z)");
    System.out.println("f(x,y,z) = (y*y*y*y*z*z*z*z+");
    System.out.println("            2*x*x*x*x*z*z*z*z+");
    System.out.println("            3*x*x*x*x*y*y*y*y+");
    System.out.println("            4+12*z*x*y+");
    System.out.println("            4*y*y*z*z*x*x)*exp(x*y*z)");
    System.out.println("");
    System.out.println("non uniform grid");
    System.out.println("");
    System.out.println("known solution, U(x,y,z), to test method");
    System.out.println("u(x,y,z) = exp(x*y*z)");
    System.out.println("xmin="+f6.format(xg[0])+
		       ", xmax="+f6.format(xg[nx-1])+
                       ", nx="+nx); 
    System.out.println("ymin="+f6.format(yg[0])+
		       ", ymax="+f6.format(yg[ny-1])+
                       ", ny="+ny); 
    System.out.println("zmin="+f6.format(zg[0])+
		       ", zmax="+f6.format(zg[nz-1])+
                       ", nz="+nz); 
    System.out.println("u(xg[0],yg[0],zg[0])="+
                       fe.format(u(xg[0],yg[0],zg[0])));
    System.out.println("f(xg[0],yg[0],zg[0])="+
                       fe.format(f(xg[0],yg[0],zg[0])));
    System.out.println("u(xg[nx-1],yg[ny-1],zg[nz-1])="+
                       fe.format(u(xg[nx-1],yg[ny-1],zg[nz-1])));
    System.out.println("f(xg[nx-1],yg[ny-1],zg[nz-1])="+
                       fe.format(f(xg[nx-1],yg[ny-1],zg[nz-1])));

    init_matrix();
    System.out.println("matrix initialized");
    System.out.println("initial matrix, left side, upper ");
    for(int i=0; i<7; i++)
    {
      for(int j=0; j<7; j++)
      {
	System.out.print(" "+fe.format(ut[i][j]));
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, right side, upper ");
    for(int i=0; i<7; i++)
    {
      for(int j=nxyz-5; j<nxyz+1; j++)
      {
	System.out.print(" "+fe.format(ut[i][j]));
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, left side, lower ");
    for(int i=nxyz-6; i<nxyz; i++)
    {
      for(int j=0; j<7; j++)
      {
	System.out.print(" "+fe.format(ut[i][j]));
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, right side, lower ");
    for(int i=nxyz-6; i<nxyz; i++)
    {
      for(int j=nxyz-5; j<nxyz+1; j++)
      { 
	  System.out.print(" "+fe.format(ut[i][j]));
      }
      System.out.println(" ");
    } 
    System.out.println(" ");
  
    new simeq(nxyz, ut, us);
    printexact();

    check();
    avgerror = error/(double)nxyz;
    System.out.println("total error="+fe.format(error)+
                       ", average error="+fe.format(avgerror)+
                       ", max error="+fe.format(maxerror));
  } // end pdenu34_eq

  int S(int i, int j, int k)
  {
      return (nz-2)*(ny-2)*(i-1) + (nz-2)*(j-1) + (k-1);
  }

    double u(double x, double y, double z) // for boundary and optional check 
  {
    return Math.exp(x*y*z);
  }

    double f(double x, double y, double z) // RHS of PDE to solve 
  {
    return (y*y*y*y*z*z*z*z+
            2.0*x*x*x*x*z*z*z*z+
            3.0*x*x*x*x*y*y*y*y+
            4.0+12.0*z*x*y+
            4.0*y*y*z*z*x*x)*Math.exp(x*y*z);
  }

  void check() // typically not known, instrumentation 
  {
    double uexact, err;
  
    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 check 

  void printexact() // typically not known 
  {
  
    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++)
        {
	  System.out.println("xg["+i+"]="+f6.format(xg[i])+
			     ", yg["+j+"]="+f6.format(yg[j])+
			     ", zg["+k+"]="+f6.format(zg[k])+
			     ", u="+fe.format(u(xg[i],yg[j],zg[k]))+
                             ", us="+fe.format(us[S(i,j,k)])+
			     ", err="+fe.format(u(xg[i],yg[j],zg[k])-
                                                us[S(i,j,k)]));
	}
      }
    }
    System.out.println(" ");
  } // end printexact 

  void init_matrix() // based on PDE
  {
    int cs, ijk, iii, jjj, kkk;
    double chk;

    cs = nxyz; // RHS, column subscript 
  
    // zero internal cells (zero based row and column)
    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;
          ut[S(i,j,k)][cs] = 0.0;
	}
    System.out.println("internal cells zeroed ");          
    for(int i=1; i<nx-1; i++) // inside boundary 
    {
      for(int j=1; j<ny-1; j++) // inside boundary 
      {
        for(int k=1; k<nz-1; k++) // inside boundary 
        {
	  ijk = S(i,j,k); // row of matrix 
          // for each term 

          // coefdxxxx * d^4 u(x,y,z)/dx^4 
          coefdxxxx = 1.0;
          new nuderiv(4,nx,i,xg,cxxxx);
          for(int ic=0; ic<nx; ic++)
	  {
            cxxxx[ic] = cxxxx[ic] * coefdxxxx;
            iii = S(ic,j,k);
	    if(ic==0 || ic==nx-1)
	      ut[ijk][cs] = ut[ijk][cs] - cxxxx[ic]*u(xg[ic],yg[j],zg[k]);
	    else ut[ijk][iii] = ut[ijk][iii] + cxxxx[ic];
	  }

          // coefdyyyy * d^4 u(x,y,z)/dy^4 
          coefdyyyy = 2.0;
          new nuderiv(4,ny,j,yg,cyyyy);
          for(int jc=0; jc<ny; jc++)
	  {
            cyyyy[jc] = cyyyy[jc] * coefdyyyy;
            iii = S(i,jc,k);
	    if(jc==0 || jc==ny-1)
	      ut[ijk][cs] = ut[ijk][cs] - cyyyy[jc]*u(xg[i],yg[jc],zg[k]);
	    else ut[ijk][iii] = ut[ijk][iii] + cyyyy[jc];
	  }

          // coefdzzzz * d^4 u(x,y,z)/dz^4 
          coefdzzzz = 3.0;
          new nuderiv(4,nz,k,zg,czzzz);
          for(int kc=0; kc<nz; kc++)
	  {
            czzzz[kc] = czzzz[kc] * coefdzzzz;
            iii = S(i,j,kc);
	    if(kc==0 || kc==nz-1)
	      ut[ijk][cs] = ut[ijk][cs] - czzzz[kc]*u(xg[i],yg[j],zg[kc]);
	    else ut[ijk][iii] = ut[ijk][iii] + czzzz[kc];
	  }

          // coefdxdydz * d^3 u(x,y)/dxdydz 
          coefdxdydz = 4.0;
          new nuderiv(1,nx,i,xg,cx);
          new nuderiv(1,ny,j,yg,cy);
          new nuderiv(1,nz,k,zg,cz);
          for(int ic=0; ic<nx; ic++)
	  {
            for(int jc=0; jc<ny; jc++)
	    {
              for(int kc=0; kc<nz; kc++)
	      {
                cxyz[ic][jc][kc] = cx[ic] * cy[jc] * cz[kc] * coefdxdydz;
                iii = S(ic,jc,kc);
	        if(ic==0 || ic==nx-1 || jc==0 || jc==ny-1 || kc==0 || kc==nz-1)
                   ut[ijk][cs] = ut[ijk][cs] - cxyz[ic][jc][kc]*
		                               u(xg[ic],yg[jc],zg[kc]);
	        else ut[ijk][iii] = ut[ijk][iii] + cxyz[ic][jc][kc];
	      }
            }
	  }

          // RHS constant 
          ut[ijk][cs] = ut[ijk][cs] + f(xg[i],yg[j],zg[k]);
          // end terms
        } // end k inside boundary
      } // end j inside boundary
    } // end i inside boundary
  } // end init_matrix 

  // check functions, not usually known
  double uxxxx(double x, double y, double z)
  {
    return y*y*y*y*z*z*z*z*Math.exp(x*y*z);
  }

  double uyyyy(double x, double y, double z)
  {
    return x*x*x*x*z*z*z*z*Math.exp(x*y*z);
  }

  double uzzzz(double x, double y,double z)
  {
    return x*x*x*x*y*y*y*y*Math.exp(x*y*z);
  }

  double uxyz(double x, double y,double z)
  {
    return Math.exp(x*y*z)+
           3.0*z*x*y*Math.exp(x*y*z)+
	   y*y*z*z*x*x*Math.exp(x*y*z);
  }

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