// corner0.java  2D uniform grid boundary value PDE 
//                known solution to test method
//              alpha = Pi/omega
//              r = sqrt(x^2+y^2)
//              theta = atan2(y,x) 0..2Pi
//              u(x,y) = r^alpha * sin(alpha*theta)         
//                                                             
//              differential equation to solve                 
//              d^2u/dx^2 + d^2u/dy^2 = 0
//                                                             
// uniform grid on rectangle -1 to 1, -1 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 corner0
{
  int nx = 11;          // includes boundary 
  int ny = 11;          // 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 = -1.0;
  double xmax =  1.0;
  double ymin = -1.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");

  corner0()
  {
    System.out.println("corner0.java running");
    System.out.println("differential equation to solve");
    System.out.println("d^2u/dx^2 + d^2u/dy^2 = 0");
    System.out.println("uniform grid on rectangle -1,1 to -1,1");
    System.out.println("known Solution, for testing method");
    System.out.println("omega = 2Pi");
    System.out.println("alpha = Pi/omega");
    System.out.println("theta = atan2(y,x) 0..2Pi");
    System.out.println("u(x,y) = r^alpha * sin(alpha*theta)");
    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)); 
    checkU( 1.0, 0.0);
    checkU( 1.0, 1.0);
    checkU( 0.0, 1.0);
    checkU(-1.0, 1.0);
    checkU(-1.0, 0.0);
    checkU(-1.0,-1.0);
    checkU( 0.0,-1.0);
    checkU( 1.0,-1.0);

    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 corner0.java");
  } // end corner0


  double u(double x, double y) // for boundary and optional check 
  {
      double omega, alpha, r, theta;

      omega = 2.0*Math.PI;
      alpha = Math.PI/omega;
      r = Math.sqrt(x*x+y*y);
      theta = Math.atan2(y,x);
      if(theta<0.0) theta = 2.0*Math.PI+theta;
      return Math.pow(r,alpha) * Math.sin(alpha*theta);
  } // end u

  void checkU(double x, double y)
  {
    double theta = Math.atan2(y,x);

    if(theta<0.0) theta = 2.0*Math.PI+theta; // 0..2Pi
    System.out.println("x="+x+", y="+y+", th="+theta+",u="+u(x,y));
  } // end checkU

  double c(double x, double y) // RHS of differential equation to solve 
  {
    return 0.0;
  } // 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("corner0_java.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing corner0_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("corner0_java.dat  not opened");
    }
    catch(Exception exception)
    {
      System.out.println("corner0_java.dat  not opened");
    }
    System.out.println("finished writing corner0_java.dat");
  }  // end write_soln2 

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