// pdenu22_eq.java  2D non uniform grid boundary value PDE  non iterative 
//                known solution to test method                         
//                 u(x,y) = x^3+2y^3+3x^2y+4xy^2+5x^2+6y^2+7x+8
//                                                                      
//                differential equation to solve                       
// d^2u/dx^2 + 2*d^2u/dy^2 + 3*d^2u/dxdy + 4*du/dx + 5*du/dy +6*u = c(x,y)
// c(x,y) = 6x^3+12y^3+18x^2y+24xy^2+57x^2+82y^2+64xy+122x+114y+110 
// non uniform grid on rectangle -1.0,-1.1 to 1.2,1.3                
//                                                                   
// set up and solve system of linear equations                       
//                                                                   
// create two dimensional array with                                 
// (nx-2)*(ny-2) rows, one for each equation                         

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

class pdenu22_eq
{
  int nx = 11;          // includes boundary 
  int ny = 10;          // includes boundary 
  double xg[] = {-1.0, -0.9, -0.8, -0.6, -0.3, 0.0,
                  0.35, 0.65, 0.85, 1.0,  1.2};
  double yg[] = {-1.1, -1.0, -0.9, -0.6, -0.2,
                  0.3,  0.8,  1.0,  1.2,  1.3};
                                 // nx-2 by ny-2  internal, non boundary cells 
  double cxx[] = new double[nx];    // nuderiv coefficients 
  double cyy[] = new double[ny];
  double cxy[][] = new double[nx][ny];
  double cx[] = new double[nx];
  double cy[] = new double[ny];
  double coefdxx, coefdyy, coefdxdy, coefdx, coefdy, coefu;

  int nxy = (nx-2)*(ny-2); // inside boundaries
  double us[] = new double[nxy];          // solution being computed 9*8 
  double ut[][] = new double[nxy][nxy+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");

  pdenu22_eq()
  {

    System.out.println("pdenu22_eq.java running");
    System.out.println("differential equation to solve");
    System.out.println("d^2u/dx^2 + 2*d^2u/dy^2 + 3*d^2u/dxdy + 4*du/dx + 5*du/dy +6*u = c(x,y)");
    System.out.println("c(x,y) = 6x^3+12y^3+18x^2y+24xy^2+57x^2+82y^2+64xy+122x+114y+110");
    System.out.println("non uniform grid on rectangle -1.0,-1.1 to 1.2,1.3");
    System.out.println("known solution, U(x,y), to test method");
    System.out.println("u(x,y) = x^3+2y^3+3x^2y+4xy^2+5x^2+6y^2+7x+8");
    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("u(xg[0],yg[0])="+f6.format(u(xg[0],yg[0])));
    System.out.println("c(xg[0],yg[0])="+f6.format(c(xg[0],yg[0])));

    init_matrix();
    System.out.println("matrix initialized");
    System.out.println("initial matrix, left side, upper ");
    for(int i=0; i<8; i++)
    {
      for(int j=0; j<8; j++)
      {
	System.out.print(" "+f6.format(ut[i][j]));
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, right side, upper ");
    for(int i=0; i<8; i++)
    {
      for(int j=65; j<73; j++)
      {
	System.out.print(" "+f6.format(ut[i][j]));
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, left side, lower ");
    for(int i=64; i<72; i++)
    {
      for(int j=0; j<8; j++)
      {
	System.out.print(" "+f6.format(ut[i][j]));
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, right side, lower ");
    for(int i=64; i<72; i++)
    {
      for(int j=65; j<73; j++)
      { 
	  System.out.print(" "+f6.format(ut[i][j]));
      }
      System.out.println(" ");
    } 
    System.out.println(" ");
  
    new simeq(nxy, ut, us);
    printexact();

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


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

  double c(double x, double y) // RHS of differential equation to solve 
  {
    return 6.0*x*x*x + 12.0*y*y*y + 18.0*x*x*y + 24.0*x*y*y + 57.0*x*x +
	   82.0*y*y + 64.0*x*y + 122.0*x + 114.0*y + 110.0;
  }

  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++)
      { 
        uexact = u(xg[i], yg[j]);
        err = Math.abs(us[(i-1)+(nx-2)*(j-1)]-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++)
      {
	  System.out.println("xg["+i+"]="+f6.format(xg[i])+
			     ", yg["+j+"]="+f6.format(yg[j])+
			     ", u="+f6.format(u(xg[i],yg[j]))+
                             ", us="+f6.format(us[(i-1)+(nx-2)*(j-1)])+
                             ", err="+fe.format(u(xg[i],yg[j])-
                                                us[(i-1)+(nx-2)*(j-1)]));
      }
    }
    System.out.println(" ");
  } // end printexact 

  void init_matrix() // based on PDE
  {
      int ik, kj, kxky, row, cs;

    cs = (nx-2)*(ny-2); // constant column subscript 
  
    // zero internal cells (zero based row and column)
    for(int i=0; i<nxy; i++)
      for(int j=0; j<nxy+1; j++)
        ut[i][j] = 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 
      {
        row = (i-1)+(nx-2)*(j-1); // row of matrix 
        // for each term 

        // coefdxx * d^2 u(x,y)/dx^2 
        coefdxx = 1.0;
        new nuderiv(2,nx,i,xg,cxx);
        if(i==1 && j==1)
	{
          System.out.println("cxx ");
	  for(int k=0; k<nx; k++)
	    System.out.println(cxx[k]);
	}
        for(int k=0; k<nx; k++)
	{
          cxx[k] = cxx[k] * coefdxx;
          kj = (k-1)+(nx-2)*(j-1);
	  if(k==0 || k==nx-1)
               ut[row][cs] = ut[row][cs] - cxx[k]*u(xg[k],yg[j]);
	  else ut[row][kj] = ut[row][kj] + cxx[k];
	}

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

        // coefdxdy * d^2 u(x,y)/dxdy 
        coefdxdy = 3.0;
        new nuderiv(1,nx,i,xg,cx);
        new nuderiv(1,ny,j,yg,cy);
        for(int kx=0; kx<nx; kx++)
	{
          for(int ky=0; ky<ny; ky++)
	  {
            cxy[kx][ky] = cx[kx] * cy[ky] * coefdxdy;
            kxky = (kx-1)+(nx-2)*(ky-1);
	    if(kx==0 || kx==nx-1 || ky==0 || ky==ny-1)
                 ut[row][cs]   = ut[row][cs]   - cxy[kx][ky]*u(xg[kx],yg[ky]);
	    else ut[row][kxky] = ut[row][kxky] + cxy[kx][ky];
          }
	}

        // coefdx * d u(x,y)/dx 
        coefdx = 4.0;
        new nuderiv(1,nx,i,xg,cx);
        for(int k=0; k<nx; k++)
	{
          cx[k] = cx[k] * coefdx;
          kj = (k-1)+(nx-2)*(j-1);
	  if(k==0 || k==nx-1)
               ut[row][cs] = ut[row][cs] - cx[k]*u(xg[k],yg[j]);
	  else ut[row][kj] = ut[row][kj] + cx[k];
	}

        // coefdy * d u(x,y)/dy 
        coefdy = 5.0;
        new nuderiv(1,ny,j,yg,cy);
        for(int k=0; k<ny; k++)
	{
          cy[k] = cy[k] * coefdy;
          ik = (i-1)+(nx-2)*(k-1);
	  if(k==0 || k==ny-1)
               ut[row][cs] = ut[row][cs] - cy[k]*u(xg[i],yg[k]);
	  else ut[row][ik] = ut[row][ik] + cy[k];
	}

        // coefu * u(x,y) 
        coefu = 6.0;
        ut[row][row] = ut[row][row] + 1.0*coefu;

        // RHS constant 
        ut[row][cs] = ut[row][cs] + c(xg[i],yg[j]);
      // end terms 
      } // end j inside boundary
    } // end i inside boundary
  } // end init_matrix 

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