// pde3b.java  3D non equal h boundary value PDE, iterative solution  
//  u(x,y,z) = x^3  + y^3  +z^3  + 2x^2y + 3xy^2 + 4z^2x + 5y + 6  
//     du/dx = 3x^2 + 0    + 0   + 4xy   + 3y^2  + 4z^2  + 0  + 0  
// d^2u/dx^2 = 6x   + 0    + 0   + 4y    + 0     + 0     + 0  + 0  
//     du/dy = 0    + 3y^2 + 0   + 2x^2  + 6xy   + 0     + 5  + 0  
// d^2u/dy^2 = 0    + 6y   + 0   + 0     + 6x    + 0     + 0  + 0  
//     du/dz = 0    + 0    +3z^2 + 0     + 0     + 8zx   + 0  + 0  
// d^2u/dz^2 = 0    + 0    +6z   + 0     + 0     + 8x    + 0  + 0  
//                                                        
// solve d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 = c(x,y,z)     
// by second order method on cube -1,-1,-1 to 1.2,1,1     

// second order numerical difference                      
// d^2u/dx^2 + d^2u/dy^2 + d^2u/dz^2 =                    
// 1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +            
// 1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z))              
// 1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) + O(h^2)     
//  = c(x,y,z) = 20x + 10y + 6z                           
// solve for new u(x,y,z) =                               
// ((u(x-hx,y,z)+u(x+hx,y,z))/hx^2 +                      
//  (u(x,y-hy,z)+u(x,y+hy,z))/hy^2 +                      
//  (u(x,y,z-hz)+u(x,y,z+hz))/hz^2 - c(x,y,z)) /          
//  (2/hx^2 + 2/hy^2 + 2/hz^2)                            


public class pde3b
{
  final int nx = 11;
  final int ny = 11;
  final int nz = 9;
  double us[][][] = new double[nx][ny][nz]; // solution being iterated        
  double ut[][][] = new double[nx][ny][nz]; // temporary for partial results 
                     // nx-2 by ny-2 by nz-2 internal, non boundary cells 
  final double xmin = -1.0;
  final double ymin = -1.0;
  final double zmin = -1.0;
  final double xmax =  1.2;
  final double ymax =  1.0;
  final double zmax =  1.0;
  double hx;   // x step = (xmax-xmin)/(nx-1) 
  double hy;   // y step = (ymax-ymin)/(ny-1) 
  double hz;   // z step = (zmax-zmin)/(nz-1) 
  double diff=1.0e100;  // sum abs difference form last iteration 
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     

  public pde3b()
  {
    int iter=0;
  
    System.out.println("pde3b.java running");
    hx = (xmax-xmin)/(nx-1); // step size in x 
    hy = (ymax-ymin)/(ny-1); // step size in y 
    hz = (zmax-zmin)/(nz-1); // step size in z 
    System.out.println("xmin="+xmin+ "xmax="+xmax+ "nx="+nx+ "hx="+hx);
    System.out.println("ymin="+ymin+ "ymax="+ymax+ "ny="+ny+ "hy="+hy);
    System.out.println("zmin="+zmin+ "zmax="+zmax+ "nz="+nz+ "hz="+hz);
  
    init_boundary();
    System.out.println("initial conditions set");
  
    while(diff>0.0001 && iter<149)
    {
      iterate();
      check();
      avgerror = error/((nx-2)*(ny-2)*(nz-2));
      System.out.println("iter="+iter+ "  diff="+diff+ 
                         "  error="+error+ "  avg="+avgerror+ "  max="+maxerror);
      iter++;
      if(iter%50==0) { printus(); printdiff(); }
    }
    printexact();
    printus();
    printdiff();
  } // end pde3b

  double u(double x, double y, double z) // solution for boundary and check 
  {
    return x*x*x + y*y*y + z*z*z + 2.0*x*x*y + 3.0*x*y*y +
           4.0*z*z*x + 5.0*y +6.0;
  }

  double c(double x, double y, double z) // from solve 
  {
    return 20.0*x + 10.0*y + 6.0*z;
  }

  double u000(int i, int j, int k) // new iterated value 
  {
    return ((us[i-1][j][k] + us[i+1][j][k])/(hx*hx) + 
            (us[i][j-1][k] + us[i][j+1][k])/(hy*hy) +
            (us[i][j][k-1] + us[i][j][k+1])/(hz*hz) -
            c(xmin+i*hx,ymin+j*hy,zmin+k*hz))       /
            (2.0/(hx*hx) + 2.0/(hy*hy) + 2.0/(hz*hz));
  }

  void init_boundary()
  {
    int i, j, k;
  
    // set boundary on six sides 
    for(i=0; i<nx; i++)
    {
      for(j=0; j<ny; j++)
      {
        us[i][j][0]    = u(xmin+i*hx, ymin+j*hy, zmin);
        us[i][j][nz-1] = u(xmin+i*hx, ymin+j*hy, zmax);
      }
    }
    for(j=0; j<ny; j++)
    {
      for(k=0; k<nz; k++)
      {
        us[0][j][k]    = u(xmin, ymin+j*hy, zmin+k*hz);
        us[nx-1][j][k] = u(xmax, ymin+j*hy, zmin+k*hz);
      }
    }
    for(i=0; i<nx; i++)
    {
      for(k=0; k<nz; k++)
      {
        us[i][0][k]    = u(xmin+i*hx, ymin, zmin+k*hz);
        us[i][ny-1][k] = u(xmin+i*hx, ymax, zmin+k*hz);
      }
    }
    // zero internal cells 
    for(i=1; i<nx-1; i++)
      for(j=1; j<ny-1; j++)
        for(k=1; k<nz-2; k++)
          us[i][j][k] = 0.0;
  } // end init_boundary 

  void iterate()
  {
    int i, j, k;
  
    diff = 0.0;
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        for(k=1; k<nz-1; k++)
        {
          ut[i][j][k] = u000(i,j,k); 
          diff = diff + Math.abs(ut[i][j][k]-us[i][j][k]);
        }
      }
    }
    // copy temp to solution 
    for(i=1; i<nx-1; i++)
      for(j=1; j<ny-1; j++)
        for(k=1; k<nz-1; k++)
          us[i][j][k] = ut[i][j][k];
  } // end iterate

  void check() // typically not known, instrumentation 
  {
    int i, j, k;
    double zexact, err;
  
    error = 0.0;
    maxerror = 0.0;
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        for(k=1; k<nz-1; k++)
        {
          zexact = u(xmin+i*hx,ymin+j*hy,zmin+k*hz);
          err = Math.abs(us[i][j][k]-zexact);
          error = error + err;
          if(err>maxerror) maxerror=err;
        }
      }
    }
  } // end check

  void printus()
  {
    int i, j, k;
  
    System.out.println("\n computed solution");
    for(i=0; i<nx; i++)
    {
      for(j=0; j<ny; j++)
      {
        for(k=0; k<nz; k++)
        {
          System.out.println("us["+i+"]["+j+"]["+k+"]= "+us[i][j][k]);
        }
      }
    }
  } // end printus

  void printexact()
  {
    int i, j, k;
  
    System.out.println("\n exact solution");
    for(i=0; i<nx; i++)
    {
      for(j=0; j<ny; j++)
      {
        for(k=0; k<nz; k++)
        {
          System.out.println("u["+i+"]["+j+"]["+k+"]="+
                             u(xmin+i*hx,ymin+j*hy,zmin+k*hz));
        }
      }
    }
  } // end printexact

  void printdiff()
  {
    int i, j, k;
  
    System.out.println("\n difference exact-computed solution");
    for(i=1; i<nx-1; i++)
    {
      for(j=1; j<ny-1; j++)
      {
        for(k=1; k<nz-1; k++)
        {
          System.out.println("diff["+i+"]["+j+"]["+k+"]= "+
                     (u(xmin+i*hx,ymin+j*hy,zmin+k*hz)-us[i][j][k]));
        }
      }
    }
  } // end printdiff
  
  public static void main (String[] args)
  {
    new pde3b();
  } // end main
} // end class pde3b of pde3b.java