// pdenu_eq.java  2D non uniform grid boundary value PDE  non iterative 
//                known solution to test method                         
//                u(x,y) = exp(x*y)*sin(x+y)                            
//                                                                      
//                differential equation to solve                       
// d^2u/dx^2 + 2*d^2u/dy^2 - 3*d^2u/dxdy + y*du/dx - x*du/dy =          
// exp(x*y)*sin(x+y)*(x*x + 2*y*y -3*x*y -3)                         
// 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 pdenu_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[11];    // nuderiv coefficients 
  double cyy[] = new double[10];
  double cxy[][] = new double[11][10];
  double cx[] = new double[11];
  double cy[] = new double[10];
  double coefdxx, coefdyy, coefdxdy, coefdx, coefdy, coefu;

  double us[] = new double[72];          // solution being computed 9*8 
  double ut[][] = new double[72][73];    // 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");


  pdenu_eq()
  {
    System.out.println("pdenu_eq.java running");
    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((nx-2)*(ny-2), ut, us);
    printexact();
    printus();
    printdiff();

    check();
    avgerror = error/((nx-2)*(ny-2));
    System.out.println("total error="+fe.format(error)+
                       ", average error="+fe.format(avgerror)+
                       ", max error="+fe.format(maxerror));
  } // end pdenu_eq


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

  double c(double x, double y) // RHS of differential equation to solve 
  {
    return Math.exp(x*y)*Math.sin(x+y)*(x*x + 2.0*y*y -3.0*x*y -3.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 printus()
  {
  
    System.out.println("computed solution");
    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(xg[j])+
			     ", us="+f6.format(us[(i-1)+(nx-2)*(j-1)]));
      }
    }
    System.out.println(" ");
  } // end printus 

  void printexact() // typically not known 
  {
  
    System.out.println("exact solution");
    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(xg[j])+
			     ", u(x,y)="+f6.format(u(xg[i],yg[j]))+
                             ", err="+fe.format(u(xg[i],yg[j])-
                                                us[(i-1)+(nx-2)*(j-1)]));
      }
    }
    System.out.println(" ");
  } // end printexact 

  void printdiff()
  {
  
    System.out.println("difference exact-computed solution");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<Math.min(ny-1,9); j++)
      {
	  System.out.print(" "+f6.format(u(xg[i],yg[j])-
                                           us[(i-1)+(nx-2)*(j-1)]));
      }
      System.out.println(" ");
    }
    System.out.println(" ");
  } // end printdiff 

  void init_matrix()
  {
      int ik, kj, kxky, ii, cs;

    cs = (nx-2)*(ny-2); // constant 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(    ii=1; ii<nx-1; ii++)
         for(int jj=1; jj<ny-1; jj++)
           ut[(i-1)+(nx-2)*(j-1)][(ii-1)+(nx-2)*(jj-1)] = 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 
      {
        ii = (i-1)+(nx-2)*(j-1); // row of matrix 
        // for each term 

        // coefdxx * d^2 u(x,y)/dx^ 
        coefdxx = 1.0;
        new nuderiv(2,nx,i,xg,cxx);
        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[ii][cs] = ut[ii][cs] - cxx[k]*u(xg[k],yg[j]);
	  else ut[ii][kj] = ut[ii][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[ii][cs] = ut[ii][cs] - cyy[k]*u(xg[i],yg[k]);
	  else ut[ii][ik] = ut[ii][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[ii][cs]   = ut[ii][cs]   - cxy[kx][ky]*u(xg[kx],yg[ky]);
	    else ut[ii][kxky] = ut[ii][kxky] + cxy[kx][ky];
          }
	}

        // coefdx * d u(x,y)/dx 
        coefdx = yg[j];
        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[ii][cs] = ut[ii][cs] - cx[k]*u(xg[k],yg[j]);
	  else ut[ii][kj] = ut[ii][kj] + cx[k];
	}

        // coefdy * d u(x,y)/dy 
        coefdy = -xg[i];
        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[ii][cs] = ut[ii][cs] - cy[k]*u(xg[i],yg[k]);
	  else ut[ii][ik] = ut[ii][ik] + cy[k];
	}

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

        // RHS constant 
        ut[ii][cs] = ut[ii][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 pdenu_eq();
  }
} // end class pdenu_eq