// pde_bihar2d_eq.java  discrete solution biharmonic 4th order
//
// solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) + 
//        2*uxx(x,y) + 2*uyy(x,y)  + 2*u(x,y) = f(x,y)
//
//  f(x,y) := 2.4 + 0.08*x*x + 0.32*y*y + 1.6*x*y + 0.02*x*x*x*x +
//            0.04*y*y*y*y - 0.08*x*x*y*y + 0.2*x*x*x*y + 0.4*y*y*y*x
//
//  in the rectangle  xmin< x <xmax and  ymin< y <ymax  ,
//
//     subject to boundary conditions of the first
//     kind, Dirichlet, ux and uy
//
//  Dirichlet boundary conditions are computed using:
//  ub(x,y) := 0.01*x*x*x*x + 0.02*y*y*y*y - 0.04*x*x*y*y +
//             0.1*x*x*x*y + 0.2*y*y*y*x - x*y +1
 
import java.io.*;
//                uses nuderiv.java and simeq.java

class pde_bihar2d_eq
{
  // i, j, nx, ii, jj, ny, xg, yg needed by discretization, available at call 
  int i, j, ii, jj, nx, ny;
  double xg[];   // X grig
  double yg[];   // Y grid
  double ug[];   // solution
  double kg[][]; // matrix
  double fg[];   // RHS
  double maxerr;

  pde_bihar2d_eq(){}; // construct the class once, run many times

  void run(double xmax, double ymax, int anx, int any)
  {
    nx = anx; // needed global
    ny = any;
    int nxy = nx*ny;
    xg = new double[nx];
    yg = new double[ny];
    double x, y, hx, hy;
    double val, err, avgerr;

    double xmin = 0.0; // problem parameters 
    double ymin = 0.0;

    kg = new double[nxy][nxy];  
    fg = new double[nxy];
    ug = new double[nxy];
    double Ua[] = new double[nxy];

    System.out.println("pde_bihar2d_eq.java running.");
    System.out.println("solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) + ");
    System.out.println("       2*uxx(x,y) + 2*uyy(x,y)  + 2*u(x,y) = f(x,y)");
    System.out.println(" ");
    System.out.println("f(x,y) := 2.4 + 0.08*x*x + 0.32*y*y + 1.6*x*y + 0.02*x*x*x*x +");
    System.out.println("          0.04*y*y*y*y - 0.08*x*x*y*y + 0.2*x*x*x*y + 0.4*y*y*y*x");
    System.out.println(" ");
    System.out.println("in the rectangle   xmin< x <xmax  and  ymin< y <ymax,");
    System.out.println(" ");
    System.out.println("   subject to boundary conditions of the first");
    System.out.println("   kind, Dirichlet.");
    System.out.println("   ub(x,y) := 0.01*x*x*x*x + 0.02*y*y*y*y - 0.04*x*x*y*y +");
    System.out.println("              0.1*x*x*x*y + 0.2*y*y*y*x - x*y + 1");
    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax);
    System.out.println("nx="+nx+", ny="+ny);
    System.out.println("x grid and analytic solution at ymin  ");
    hx = (xmax-xmin)/(double)(nx-1);
    for(i=0; i<nx; i++)
    {
      xg[i] = xmin+i*hx;
      Ua[i] = ub(xg[i],ymin); // temp 
      System.out.println("i="+i+", Ua("+xg[i]+")="+Ua[i]);
    }  
    System.out.println(" ");
    System.out.println("y grid and analytic solution at xmin ");
    hy = (ymax-ymin)/(double)(ny-1);
    for(ii=0; ii<ny; ii++)
    {
      yg[ii] = ymin+ii*hy;
      Ua[ii] = ub(xmin, yg[ii]); // temp 
      System.out.println("ii="+ii+", Ua("+yg[ii]+")="+Ua[ii]);
    }  
    System.out.println(" ");
    // put solution in Ua vector 
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      for(ii=0; ii<ny; ii++)
      {
	y = yg[ii];
        Ua[i*ny+ii] = ub(x,y);
        // System.out.println("solution at i="+i+",x="+x+", ii="+ii+",y="+y+
        //                    " is"+Ua[i*ny+ii]);
      }
    }
    System.out.println(" ");

    // initialize k and f 
    for(i=0; i<nx; i++)
    {
      for(j=0; j<nx; j++)
      {
        for(ii=0; ii<ny; ii++)
        {
          for(jj=0; jj<ny; jj++)
          {
            kg[(i*ny+ii)][(j*ny+jj)] = 0.0;
	  }
        }
      }
    }
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        fg[i*ny+ii] = 0.0;
      }
    }
    // set boundary conditions, four sides 
    for(i=0; i<nx; i++)
    {
      x = xg[i];
      ii=0;
      y = yg[ii];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
      // System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
      //                    " is "+fg[i*ny+ii]);
      ii=ny-1;
      y = yg[ii];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
      // System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
      //                    " is "+fg[i*ny+ii]);
    }
    for(ii=0; ii<ny; ii++)
    {
      y = yg[ii];
      i=0;
      x = xg[i];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
      // System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
      //                    " is "+fg[i*ny+ii]);
      i=nx-1;
      x = xg[i];
      kg[(i*ny+ii)][(i*ny+ii)] = 1.0;
      fg[i*ny+ii] = ub(x,y);
      // System.out.println("boundary i="+i+",x="+x+", ii="+ii+",y="+y+
      //                    " is "+fg[i*ny+ii]);
    }

    // compute value for fg[], RHS  
    for(i=1; i<nx-1; i++)
    {
      x = xg[i];
      for(ii=1; ii<ny-1; ii++)
      {
	y=yg[ii];
	fg[i*ny+ii] = F(x,y);
      }
    }

    // compute entries in kg (stiffness matrix) discretization matrix 
    System.out.println("compute discretization matrix  ");
    for(i=1; i<nx-1; i++) // not boundary rows, done above
    {
      for(ii=1; ii<ny-1; ii++) // not boundary rows, done above
      {
        // compute discrete derivatives around (i,ii) DOF  
	discretize(); // everything global
                      // sets kg[(i*ny+ii)][(j*ny+jj)...] 
      } // end ii loop 
    } // end i loop 

    System.out.println("k computed discrete matrix");
    System.out.println("f computed RHS");

    new simeq(kg, fg, ug);

    write_soln(); // to file x y soln

    check_soln(xmax, ymax, nx, ny); // independent of nowing analytic solution

    avgerr = 0.0;
    maxerr = 0.0;
    System.out.println("ug computed, Ua analytic, error  ");
    for(i=0; i<nx; i++)
    {
      for(ii=0; ii<ny; ii++)
      {
        err = Math.abs(ug[i*ny+ii]-Ua[i*ny+ii]);
        if(err>maxerr) maxerr = err;
        avgerr = avgerr + err;
        System.out.println("ug["+i+","+ii+"]="+ug[i*ny+ii]+", Ua="+Ua[i*ny+ii]+
                           ", err="+(ug[i*ny+ii]-Ua[i*ny+ii]));
      }
    }
    System.out.println("xmax="+xmax+", ymax="+ymax+", nx="+nx+", ny="+ny);
    System.out.println("                  maxerr="+maxerr+", avgerr="+
                       avgerr/(double)(nx*ny));
    System.out.println("\n");
  } // end pde_bihar2d_eq constructor

  // PDE problem definition functions  f and ub

  double F(double x, double y)
  {
    return 2.4 + 0.08*x*x + 0.32*y*y + 1.6*x*y + 0.02*x*x*x*x +
           0.04*y*y*y*y - 0.08*x*x*y*y + 0.2*x*x*x*y + 0.4*y*y*y*x;
  }

  double ub(double x, double y)
  {
    return 0.01*x*x*x*x + 0.02*y*y*y*y - 0.04*x*x*y*y +
	   0.1*x*x*x*y + 0.2*y*y*y*x - x*y + 1.0;
  }

  // solve  uxxxx(x,y) + 2*uxxyy(x,y) + uyyyy(x,y) + 
  //        2*uxx(x,y) + 2*uyy(x,y)  + 2*u(x,y) = f(x,y)
  void discretize()  // everything global i,ii (x,y) DOF
                     // sets kg[(i*ny+ii)][(j*ny+jj)...]
  {
      if(i<1 || i>nx-1 || ii<1 || ii>ny-1) // primeter is boundary
    {
      System.out.println("discretize problem i="+i+", ii="+ii);
      return;
    }
    if(i==1)
    {
      if(ii==1)         dset(1, 1); // offset is negative of x,y point
      else if(ii==ny-2) dset(1, 3);
      else              dset(1, 2);
    }
    else if(ii==1)
    {
      if(i==1)          dset(1, 1);
      else if(i==nx-2)  dset(3, 1);
      else              dset(2, 1);
    }
    else if(i==nx-2)
    {
      if(ii==1)         dset(3, 1);
      else if(ii==ny-2) dset(3, 3);
      else              dset(3, 2);
    }
    else if(ii==ny-2)
    {
      if(i==1)          dset(1, 3);
      else if(i==nx-2)  dset(3, 3);
      else              dset(2, 3);
    }
    else                dset(2, 2); // center case
  } // end discretize

  void dset(int k, int kk) // x point, y point negative for offset
  {
    double cxxxx[] = new double[5];
    double cyyyy[] = new double[5];
    double cxx  [] = new double[5];
    double cyy  [] = new double[5];
    double cxxyy[][] = new double[5][5];
    double pxg  [] = new double[5];
    double pyg  [] = new double[5];

    if(i-k<0 || i-k+4>nx-1 || ii-kk<0 || ii-kk+4>ny-1)
    {
      System.out.println("dset problem k="+k+", kk="+kk+
                         "  for  i="+i+", ii="+ii);
      return;
    }
    System.out.println("dset k="+k+", kk="+kk+", row="+(i*ny+ii));

    for(int m=0; m<5; m++)
    {
      pxg[m] = xg[i-k+m];  // offset
      pyg[m] = yg[ii-kk+m];
    }
    new nuderiv(4, 5, k, pxg, cxxxx);
    new nuderiv(2, 5, k, pxg, cxx);
    new nuderiv(4, 5, kk, pyg, cyyyy);
    new nuderiv(2, 5, kk, pyg, cyy);
    for(int mx=0; mx<5; mx++)
    {
      for(int my=0; my<5; my++) cxxyy[mx][my] = cxx[mx]*cyy[my];
    } 

    for(int mx=0; mx<5; mx++) // process all 25 cells
    {
      for(int my=0; my<5; my++)
      {
      // optional, if boundary not in A and U
      //if(my==ii && (i-k+mx==0 || i-k+mx==nx-1)) // x boundary
      //{
      //  fg[i*ny+ii] -= ub(xg[i-k+mx],yg[ii-kk+my])*cxxxx[mx];
      //  fg[i*ny+ii] -= 2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cxx[mx];
      //  if(i==1 && ii==1)
      //  {
      //    System.out.println("bound i="+i+",ii="+ii+", mx="+mx+", my="+my);
      //    System.out.println("RHS cxxxx- "+
      //		       (ub(xg[i-k+mx],yg[ii-kk+my])*cxxxx[mx]));
      //    System.out.println("RHS cxx- "+
      //                       (2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cxx[mx]));
      //  }
      //}
      //else if(mx==i && (ii-kk+my==0 || ii-kk+my==ny-1)) // y boundary
      //{
      //  fg[i*ny+ii] -= ub(xg[i-k+mx],yg[ii-kk+my])*cyyyy[my];
      //  fg[i*ny+ii] -= 2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cyy[my];
      //  if(i==1 && ii==1)
      //  {
      //    System.out.println("bound i="+i+",ii="+ii+", mx="+mx+", my="+my);
      //    System.out.println("RHS cyyyy- "+
      //                       (ub(xg[i-k+mx],yg[ii-kk+my])*cyyyy[my]));
      //    System.out.println("RHS cyy- "+
      //                       (2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cyy[my]));
      //  }
      //}
      //else if(i-k+mx==0 || i-k+mx==nx-1 || ii-kk+my==0 || ii-kk+my==ny-1)
      //{                                                     // xxyy boundary
      //  fg[i*ny+ii] -= 2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cxxyy[mx][my];
      //  if(i==1 && ii==1)
      //  {
      //    System.out.println("bound i="+i+",ii="+ii+", mx="+mx+", my="+my);
      //    System.out.println("RHS cxxyy- "+
      //                      (2.0*ub(xg[i-k+mx],yg[ii-kk+my])*cxxyy[mx][my]));
      //  }
      //}
      //else
      //{
	  if(kk==my) // Uxxxx and Uxx
	  {
	    kg[(i*ny+ii)][(i-k+mx)*ny+ii] += cxxxx[mx];
	    kg[(i*ny+ii)][(i-k+mx)*ny+ii] += 2.0*cxx[mx];
	  }
	  if(k==mx) // Uyyyy and Uyy
	  {
	    kg[(i*ny+ii)][i*ny+(ii-kk+my)] += cyyyy[my];
	    kg[(i*ny+ii)][i*ny+(ii-kk+my)] += 2.0*cyy[my];
	  }
          // Uxxyy
	  kg[(i*ny+ii)][(i-k+mx)*ny+(ii-kk+my)] += 2.0*cxxyy[mx][my];
      //}
      } // end my
    } // end mx
    kg[(i*ny+ii)][(i*ny+ii)] += 2.0; // U
  } // end dset

  void write_soln()
  {
    try
    {
      FileWriter fout = new FileWriter("pde_bihar2d_eq.dat");
      BufferedWriter fileout = new BufferedWriter(fout);
      System.out.println("writing pde_bihar2d_eq.dat");

      for(int i=0; i<nx; i++)
      {
        for(int ii=0; ii<ny; ii++)
        {
          fileout.write(" "+xg[i]+" "+yg[ii]+" "+ug[i*ny+ii]+" \n");
        } // end y
        fileout.write(" \n");
      } // end x
      fileout.close();
      System.out.println("finished writing pde_bihar2d_eq.dat");
    }
    catch(FileNotFoundException exception)
    {
      System.out.println("pde_bihar2d_eq.dat  not opened");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar2d_eq.java output-file error");
    } // end try
  } // end write_soln

  void check_soln(double xmax, double ymax, int nx, int ny)
                  // for all interior points
  {

    double x, y, err, maxerr, avgerr, rmserr;
    double sumerr, sumsqerr, eval_deriv;
    double cxx[][] = new double[nx][nx];
    double cxxxx[][] = new double[nx][nx];
    double cyy[][] = new double[ny][ny];
    double cyyyy[][] = new double[ny][ny];
    double tcxxyy[] = new double[ny]; // for partial uxxyy
    int cnt = 0;

    for(int i=1; i<nx-1; i++)
    {
      new nuderiv(2,nx,i,xg,cxx[i]);
      new nuderiv(4,nx,i,xg,cxxxx[i]);
    }
    for(int ii=1; ii<ny-1; ii++)
    {
      new nuderiv(2,ny,ii,yg,cyy[ii]);
      new nuderiv(4,ny,ii,yg,cyyyy[ii]);
    }

    sumerr = 0.0;
    sumsqerr = 0.0;
    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // full equation
    {
      x = xg[i];
      for(int j=1; j<ny-1; j++)
      {
	y = yg[j];
	eval_deriv = 2.0*ug[i*ny+j]; // PDE u term
        
        for(int k=0; k<nx; k++)
        {
	  eval_deriv += 2.0*cxx[i][k]*ug[k*ny+j]; // PDE  uxx term
	  eval_deriv += cxxxx[i][k]*ug[k*ny+j];   // PDE  uxxxx term
	} 
        for(int k=0; k<ny; k++)
        {
	  eval_deriv += 2.0*cyy[j][k]*ug[i*ny+k];  // PDE uyy term
	  eval_deriv += cyyyy[j][k]*ug[i*ny+k];    // PDE uyyyy term
	}
        //  + 2*uxxyy
        for(int jj=0; jj<ny; jj++) // walk jj in y, computing uxx
	{
	  tcxxyy[jj] = 0.0; // numeric  uxx part of uxxyy
          for(int k=0; k<nx; k++)
          {
	    tcxxyy[jj] += cxx[i][k]*ug[k*ny+jj];
	  }
	}
        for(int k=0; k<ny; k++)  // uyy part of uxxyy
        {
	  eval_deriv += 2.0*cyy[j][k]*tcxxyy[k];  // PDE uxxyy term
	}

        err = eval_deriv-F(x,y);
        sumerr += Math.abs(err);
        sumsqerr += err*err;
        maxerr = Math.max(Math.abs(err),maxerr);
        cnt++;  
      } // end y
    } // end x
    rmserr = Math.sqrt(sumsqerr/(double)cnt);
    avgerr = sumerr/(double)cnt;
    System.out.println("check_soln values");
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);
  } // end check_soln 

  public static void main (String[] args)
  {
    pde_bihar2d_eq F2 = new pde_bihar2d_eq();
  
    F2.run(4.0, 4.0, 9, 9); // xmax, ymax, nx, ny

    try
    {
      Runtime runtime = Runtime.getRuntime();
      Process proc = runtime.exec("gnuplot -persist pde_bihar2d_eq.plot");
    }
    catch(Exception exception) // catches all exceptions
    {
      System.out.println("pde_bihar2d_eq.java runtime error");
    } // end try
  }
} // end class pde_bihar2d_eq