// pde22_eq.java  2D uniform grid boundary value PDE
//                known solution to test method
//     u(x,y) = 2^16 *x^4 * (1-x)^4 * y^4 * (1-y)^4;  zero on boundary 
//                                           hump, max 1.0 at  (0.5,0.5)
//     differential equation to solve
//     d^2u/dx^2 + d^2u/dy^2 = c(x,y)
//     Uxx(x,y) = 2^16 * (
//                        12*x^2 * (1-x)^4 * y^4 * (1-y)^4 - 
//                        32*x^3 * (1-x)^3 * y^4 * (1-y)^4 + 
//                        12*x^4 * (1-x)^2 * y^4 * (1-y)^4)  
//     Uyy(x,y) = 2^16 * (
//                        12*x^4 * (1-x)^4 * y^2 * (1-y)^4 - 
//                        32*x^4 * (1-x)^4 * y^3 * (1-y)^3 + 
//                        12*x^4 * (1-x)^4 * y^4 * (1-y)^2)  
//     c(x,y) = Uxx(x,y) + Uyy(x,y);
// 
// uniform grid on rectangle 0 to 1, 0 to 1                    
// set up and solve system of (nx-2)*(ny-2) linear equations                 
//
// uses nuderiv.java
// uses simeq.java
import java.text.*;
import java.io.*;

class pde22_eq
{
  int nx = 10;          // includes boundary 
  int ny = 10;          // includes boundary 
  double xg[] = new double[nx];
  double yg[] = new double[ny];
  double cxx[] = new double[nx];    // nuderiv coefficients 
  double cyy[] = new double[ny];
                                 // nx-2 by ny-2  internal, non boundary cells 
  int neqn = (nx-2)*(ny-2);
  int cs = neqn;
  double us[] = new double[neqn];               // solution being computed 
  double ut[][] = new double[neqn][neqn+1];     // equations incl RHS 
  double xmin = 0.0;
  double xmax = 1.0;
  double ymin = 0.0;
  double ymax = 1.0;
  double hx, hy;
  double error;         // sum of absolute exact-computed         
  double avgerror;      // average error                          
  double maxerror;      // maximum cell error                     
  double time_start;
  double time_now;
  double time_last;
  DecimalFormat fe = new DecimalFormat("0.00E00");
  DecimalFormat f6 = new DecimalFormat("0.00000");

  pde22_eq()
  {
    System.out.println("pde22_eq.java running");
    System.out.println("differential equation to solve");
    System.out.println("d^2u/dx^2 + d^2u/dy^2 = c(x,y)");
    System.out.println("Uxx(x,y) = 2^16 * (");
    System.out.println("            12*x^2 * (1-x)^4 * y^4 * (1-y)^4 -");
    System.out.println("            32*x^3 * (1-x)^3 * y^4 * (1-y)^4 +");
    System.out.println("            12*x^4 * (1-x)^2 * y^4 * (1-y)^4 )");
    System.out.println("Uyy(x,y) = 2^16 * (");
    System.out.println("            12*x^4 * (1-x)^4 * y^2 * (1-y)^4 -");
    System.out.println("            32*x^4 * (1-x)^4 * y^3 * (1-y)^3 +");
    System.out.println("            12*x^4 * (1-x)^4 * y^4 * (1-y)^2 )");
    System.out.println("c(x,y) = Uxx(x,y) + Uyy(x,y);");
    System.out.println("uniform grid on rectangle 0,1 to 0,1");
    System.out.println("known Solution, for testing method");
    System.out.println("u(x,y) = 2^16*x^4*(1-x)^4*y^4*(1-y)4");
    System.out.println("zero on boundary, hump 1.0 at x=0.5, y=0.5 ");
    System.out.println(" ");

    time_start = (double)System.currentTimeMillis()/1000.0;
    time_last = time_start;
    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);

    for(int i=0; i<nx; i++)
    {
      xg[i] = xmin+(double)i*hx;
      System.out.println("xg("+i+")="+xg[i]);
    }
    for(int j=0; j<ny; j++)
    {
      yg[j] = ymin+(double)j*hy;
      System.out.println("yg("+j+")="+yg[j]);
    }

    System.out.println("xmin="+f6.format(xg[0])+
		       ", xmax="+f6.format(xg[nx-1])+
                       ", nx="+nx+", hx="+f6.format(hx)); 
    System.out.println("ymin="+f6.format(yg[0])+
		       ", ymax="+f6.format(yg[ny-1])+
                       ", ny="+ny+", hy="+f6.format(hy)); 
    System.out.println("u(0.5,0.5)="+
                       f6.format(u(0.5,0.5)));
    System.out.println("c(0,5,0.5)="+
                       f6.format(c(0.5,0.5)));

    init_matrix();
    System.out.println("matrix initialized");
    System.out.println("initial matrix, left side, upper ");
    // printB();
    check_rows(); 
    new simeq(neqn, ut, us);
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("time to compute solution = "+(time_now-time_last)+
                       " seconds ");
    time_last = time_now;
    check_soln();
    print();

    check();
    avgerror = error/(double)neqn;
    System.out.println("total error="+fe.format(error)+
                       ", average error="+fe.format(avgerror)+
                       ", max error="+fe.format(maxerror));
    write_soln2(us);
    time_now = (double)System.currentTimeMillis()/1000.0;
    System.out.println("total time = "+(time_now-time_start)+
                       " seconds ");
    System.out.println("finished pde22_eq.java");
  } // end pde22_eq


  double u(double x, double y) // for boundary and optional check 
  {
    return 65536.0*  // 2^4*p  p=4 
           x*x*x*x*(1.0-x)*(1.0-x)*(1.0-x)*(1.0-x)*
           y*y*y*y*(1.0-y)*(1.0-y)*(1.0-y)*(1.0-y);
  } // end u

  double c(double x, double y) // RHS of differential equation to solve 
  {
    return 65536.0*(
           12.0*x*x*    (1.0-x)*(1.0-x)*(1.0-x)*(1.0-x)*y*y*y*y*(1.0-y)*(1.0-y)*(1.0-y)*(1.0-y)-
           32.0*x*x*x*  (1.0-x)*(1.0-x)*(1.0-x)*        y*y*y*y*(1.0-y)*(1.0-y)*(1.0-y)*(1.0-y)+
           12.0*x*x*x*x*(1.0-x)*(1.0-x)*                y*y*y*y*(1.0-y)*(1.0-y)*(1.0-y)*(1.0-y)+
           12.0*x*x*x*x*(1.0-x)*(1.0-x)*(1.0-x)*(1.0-x)*y*y*    (1.0-y)*(1.0-y)*(1.0-y)*(1.0-y)-
           32.0*x*x*x*x*(1.0-x)*(1.0-x)*(1.0-x)*(1.0-x)*y*y*y*  (1.0-y)*(1.0-y)*(1.0-y)+
           12.0*x*x*x*x*(1.0-x)*(1.0-x)*(1.0-x)*(1.0-x)*y*y*y*y*(1.0-y)*(1.0-y));
  } // end c

  int S(int i, int j) // zero based index includes boundary 
  {                           // this is row of ut[] 
    return (i-1)*(ny-2)+(j-1); // index into us and ut, RHS, no boundary 
  } // end S 

  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[S(i,j)]-uexact);
        error = error + err;
        if(err>maxerror) maxerror=err;
      }
    }
  } // end check 

  void check_rows() // only for(checking 
  {	// row is in ut solution coordinates 
    int row, col;
    double sum;

    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
        row = S(i,j);
        sum = 0.0;
        for(int ii=1; ii<nx-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
	  {
            col = S(ii,jj);
            sum = sum + ut[row][col]*u(xg[ii],yg[jj]); // u is check solution 
	  } // ii 
        } // jj 
        sum = sum - ut[row][cs];
        if(Math.abs(sum)>0.001)
        {
          System.out.println("row i="+i+", j="+j+", is bad err="+sum);
        }
      } // j 
    } //i 
  } // end check_rows 

  void print() // typically not known 
  {
    double error = 0.0;
    double max_error = 0.0;
    double avg_error = 0.0;
  
    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++)
      {
	error = u(xg[i],yg[j]) - us[S(i,j)];
        avg_error = avg_error + Math.abs(error);
        if(Math.abs(error)>max_error) max_error = Math.abs(error);

	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[S(i,j)])+
                           ", err="+fe.format(error));
      }
    }
    System.out.println("max_error="+max_error+", avg_error="+
                       avg_error/(double)neqn);
    System.out.println(" ");
  } // end print 

  void printB()
  {
    System.out.println("matrix B includes RHS");
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int ii=1; ii<nx-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
	  {
            System.out.println("i="+i+", j"+j+", ii="+ii+", jj="+jj+", ut="+
		               ut[S(i,j)][S(ii,jj)]);
	  } // jj 
        } // ii 
        System.out.println("RHS=                         "+ut[S(i,j)][cs]);
      } // j 
    } // i 
    System.out.println("");
    System.out.println("");
  } //  end printB 

  void init_matrix()
  {
    int row, kj, ik;

    // cs is RHS, constant column subscript 
    // zero internal cells 
    for(int i=1; i<nx-1; i++)
    {
      for(int j=1; j<ny-1; j++)
      {
        for(int ii=1; ii<nx-1; ii++)
        {
          for(int jj=1; jj<ny-1; jj++)
	  {
	      ut[S(i,j)][S(ii,jj)] = 0.0;
          }
        }
        ut[S(i,j)][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++)
      {
        row = S(i,j); // equation 
        // for(each term 

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

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

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

  void check_soln()  // check when solution unknown 
  {
    double smaxerr, err, x, y, uxx, uyy;

    //  ux(x,y) + uy(x,y) = c(x,y) 
    System.out.println("check_soln against PDE");
    smaxerr = 0.0;

    for(int i=1; i<nx-1; i++)  // through dof equations 
    {
      new nuderiv(2,nx,i,xg,cxx);
      x = xg[i];
      for(int j=1; j<ny-1; j++)
      {
        new nuderiv(2,ny,j,yg,cyy);
        y = yg[j];
        err = -c(x,y); // add other side of equation terms 
        uxx = 0.0;
        for(int k=0; k<nx; k++) // compute numerical derivative at (i,j) 
        {
          if(k==0 || k==nx-1)
	  {
            uxx = uxx + cxx[k]*u(xg[k],yg[j]); // boundary 
	  }
          else
	  {
            uxx = uxx + cxx[k]*us[S(k,j)];
         }
       }
        uyy = 0.0;
        for(int k=0; k<ny; k++)
        {
          if(k==0 || k==ny-1)
	  {
            uyy = uyy + cyy[k]*u(xg[i],yg[k]); // boundary 
	  }
          else
	  {
            uyy = uyy + cyy[k]*us[S(i,k)];
	  }
        }
        err = err + uxx + uyy;
        if(Math.abs(err) > smaxerr) smaxerr = Math.abs(err);
      } // ii y 
    } // i x 

    System.out.println("check soln against PDE max error="+smaxerr);
  } // end Check_Soln 

  // write_soln2  for gnuplot 
  void  write_soln2(double Ux[])
  {
    try
    {
      FileWriter fout = new FileWriter("pde22_eq_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde22_eq_java.dat");

      for(int i=0; i<nx; i++)
      {
        for(int j=0; j<ny; j++)
        {
          if(i==0 || i==nx-1 || j==0 || j==ny-1)
          {
            fileout.write(xg[i]+" "+yg[j]+" "+u(xg[i],yg[j])+"\n"); // boundary 
          }
          else
          { 
            fileout.write(xg[i]+" "+yg[j]+" "+Ux[S(i,j)]+"\n"); // computed solution 
          }
        } // j
        fileout.write("\n"); // blank line 
      } // i 
      fileout.close();
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde22_eq_java.dat  not opened");
    }
    catch(Exception exception)
    {
      System.out.println("pde22_eq_java.dat  not opened");
    }
    System.out.println("finished writing pde22_eq_java.dat");
  }  // end write_soln2 

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