// pdeCnu22_eq.java  2D non uniform grid boundary value PDE  non iterative 
//                   needs Complex.java Cnuderiv.java Csimeq.java Cinvert.java
//                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 complex linear equations
//                                                                   
// create two dimensional array with                                 
// (nx-2)*(ny-2) rows, one for each equation                         

// uses Complex.java
// uses Cnuderiv.java
// uses Cinverse.java
// uses Csimeq.java

class pdeCnu22_eq
{
  Complex Czero = new Complex(0.0,0.0);
  Complex Cone  = new Complex(1.0,0.0);

  int nx = 11;          // includes boundary 
  int ny = 10;          // includes boundary 

  Complex xg[] = {new Complex(-1.0,-0.9), 
                  new Complex(-0.9,-0.8),
                  new Complex(-0.8,-0.7),
                  new Complex(-0.6,-0.4),
                  new Complex(-0.3,-0.2),
                  new Complex( 0.0,0.0),
                  new Complex( 0.35,0.4),
                  new Complex( 0.65,0.7),
                  new Complex( 0.85,0.9),
		  new Complex( 1.0,1.2),
		  new Complex( 1.2,1.4)};
  Complex yg[] = {new Complex(-1.1,-0.8),
                  new Complex(-1.0,-0.7),
                  new Complex(-0.9,-0.6),
                  new Complex(-0.6,-0.4),
                  new Complex(-0.2,-0.1),
                  new Complex( 0.3,0.5),
                  new Complex( 0.8,0.7),
                  new Complex( 1.0,1.0),
                  new Complex( 1.2,1.4),
		  new Complex( 1.3,1.5)};
                                 // nx-2 by ny-2  internal, non boundary cells 
  Complex cxx[] = new Complex[nx];    // nuderiv coefficients 
  Complex cyy[] = new Complex[ny];
  Complex cxy[][] = new Complex[nx][ny];
  Complex cx[] = new Complex[nx];
  Complex cy[] = new Complex[ny];
  Complex coefdxx = new Complex(1.0,0.0); // from problem definition
  Complex coefdyy = new Complex(2.0,0.0);
  Complex coefdxdy = new Complex(3.0,0.0);
  Complex coefdx = new Complex(4.0,0.0);
  Complex coefdy = new Complex(5.0,0.0);
  Complex coefu = new Complex(6.0,0.0);

  int nxy = (nx-2)*(ny-2); // inside boundary
  Complex us[] = new Complex[nxy];  // solution being computed 9*8 
  Complex ut[][] = new Complex[nxy][nxy+1];    // temporary equations 
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     

  pdeCnu22_eq()
  {

    System.out.println("pdeCnu22_eq.java running");
    System.out.println("complex 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 complex rectangle (-1.0.-0.9),(-1.1,0.8)");
    System.out.println(" to (1.2,1.4),(1.3,1.5)");
    System.out.println("known complex value solution, U(x,y), x,y complex,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="+xg[0]+", xmax="+xg[nx-1]+", nx="+nx); 
    System.out.println("ymin="+yg[0]+", ymax="+yg[ny-1]+", ny="+ny); 
    System.out.println("U(xg[0],yg[0])="+u(xg[0],yg[0]));
    System.out.println("c(xg[0],yg[0])="+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<4; i++)
    {
      for(int j=0; j<4; j++)
      {
	System.out.print(" "+ut[i][j]);
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, right side, upper ");
    for(int i=0; i<1; i++)
    {
      for(int j=72; j<73; j++)
      {
	System.out.print(" "+ut[i][j]);
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, left side, lower ");
    for(int i=71; i<72; i++)
    {
      for(int j=0; j<1; j++)
      {
	System.out.print(" "+ut[i][j]);
      }
      System.out.println(" ");
    } 
    System.out.println("initial matrix, right side, lower ");
    for(int i=71; i<72; i++)
    {
      for(int j=72; j<73; j++)
      { 
	  System.out.print(" "+ut[i][j]);
      }
      System.out.println(" ");
    } 
    System.out.println(" ");
  
    new Csimeq(nxy, ut, us);
    System.out.println("Complex system of simultaneous equations solved ");
    System.out.println(" ");
    printexact();

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


  Complex u(Complex x, Complex y) // for boundary and optional check 
  {
    Complex t1 = new Complex(1.0,0.0);
    Complex t2 = new Complex(2.0,0.0);
    Complex t3 = new Complex(3.0,0.0);
    Complex t4 = new Complex(4.0,0.0);
    Complex t5 = new Complex(5.0,0.0);
    Complex t6 = new Complex(6.0,0.0);
    Complex t7 = new Complex(7.0,0.0);
    Complex t8 = new Complex(8.0,0.0);
    t1 = x.multiply(x.multiply(x));
    t2 = t2.multiply(y.multiply(y.multiply(y))); 
    t3 = t3.multiply(x.multiply(x.multiply(y)));
    t4 = t4.multiply(x.multiply(y.multiply(y)));
    t5 = t5.multiply(x.multiply(x));
    t6 = t6.multiply(y.multiply(y));
    t7 = t7.multiply(x);
    return t1.add(t2.add(t3.add(t4.add(t5.add(t6.add(t7.add(t8)))))));
  }

  Complex c(Complex x, Complex y) // RHS of differential equation to solve 
  {
    Complex t1 = new Complex( 6.0,0.0);
    Complex t2 = new Complex(12.0,0.0);
    Complex t3 = new Complex(18.0,0.0);
    Complex t4 = new Complex(24.0,0.0);
    Complex t5 = new Complex(57.0,0.0);
    Complex t6 = new Complex(82.0,0.0);
    Complex t7 = new Complex(64.0,0.0);
    Complex t8 = new Complex(122.0,0.0);
    Complex t9 = new Complex(114.0,0.0);
    Complex t10 = new Complex(110.0,0.0);
    t1 = t1.multiply(x.multiply(x.multiply(x)));
    t2 = t2.multiply(y.multiply(y.multiply(y))); 
    t3 = t3.multiply(x.multiply(x.multiply(y)));
    t4 = t4.multiply(x.multiply(y.multiply(y)));
    t5 = t5.multiply(x.multiply(x));
    t6 = t6.multiply(y.multiply(y));
    t7 = t7.multiply(x.multiply(y));
    t8 = t8.multiply(x);
    t9 = t9.multiply(y);
    return t1.add(t2.add(t3.add(t4.add(t5.add(t6.add(t7.add(t8.add(t9.add(t10)))))))));
    //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 
  {
    Complex uexact;
    double 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 = (us[(i-1)+(nx-2)*(j-1)].subtract(uexact)).abs();
        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+"]="+xg[i]+
			     ", yg["+j+"]="+yg[j]+
			     ", u="+u(xg[i],yg[j])+
                             ", us="+us[(i-1)+(nx-2)*(j-1)]+
                             ", err="+(u(xg[i],yg[j]).subtract(us[(i-1)+(nx-2)*(j-1)]).abs()));
      }
    }
    System.out.println(" ");
  } // end printexact 

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

    cs = nxy; // 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] = new Complex(0.0,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
        // System.out.println("i="+i+", j="+j+", row ="+row); 
        // for each term 

        // coefdxx * d^2 u(x,y)/dx^2 
        new Cnuderiv(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].multiply(coefdxx);
          kj = (k-1)+(nx-2)*(j-1);  // ?? k-1
	  if(k==0 || k==nx-1)
               ut[row][cs] = ut[row][cs].subtract(cxx[k].multiply(u(xg[k],yg[j])));
	  else ut[row][kj] = ut[row][kj].add(cxx[k]);
	}

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

        // coefdxdy * d^2 u(x,y)/dxdy 
        new Cnuderiv(1,nx,i,xg,cx);
        new Cnuderiv(1,ny,j,yg,cy);
        for(int kx=0; kx<nx; kx++)
	{
          for(int ky=0; ky<ny; ky++)
	  {
	    cxy[kx][ky] = cx[kx].multiply(cy[ky].multiply(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].subtract(cxy[kx][ky].multiply(u(xg[kx],yg[ky])));
	    else ut[row][kxky] = ut[row][kxky].add(cxy[kx][ky]);
          }
	}

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

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

        // coefu * u(x,y) 
        ut[row][row] = ut[row][row].add(coefu);

        // RHS constant 
        ut[row][cs] = ut[row][cs].add(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 pdeCnu22_eq();
    System.out.println(" ");
    System.out.println("pdeCnu22_eq.java  finished");
  }
} // end class pdeCnu22_eq