// test_bihar2d.java  test setup for a specific Biharmonic PDE
//                    much of the code below is just to do plotting
//                    auxiliary tests for derivative coefficients
//
// 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, the exterior normal derivatives
//     of u must be specified at the boundary.
//
//  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
//
//  Thus the derivatives in x and y are:
//
//  ux(x,y) := 0.04*x*x*x - 0.08*x*y*y + 0.3 *x*x*y + 0.2*y*y*y - y
//
//  uy(x,y) := 0.08*y*y*y - 0.8*x*x*y + 0.1*x*x*x + 0.6*x*y*y -x
//
//  continuing, for testing, thus additional derivatives are:
//  uxx(x,y) := 0.12*x*x -0.08*y*y + 0.6*x*y
//  uxxx(x,y) := 0.24*x + 0.6*y
//  uxxxx(x,y) := 0.24
//  uyy(x,y) := 0.24*y*y -0.08*x*x + 1.2*x*y
//  uyyy(x,y) := 0.48*y + 1.2*x
//  uyyyy(x,y) := 0.48
//  uxxyy(x,y) := -0.16
//

import java.io.*;
import java.*;
import java.awt.*;
import java.awt.event.*;
//                       needs nderiv.java 

public class test_bihar2d extends Frame
{
  // mouse and display variables
  int xm, ym, bm; // mouse position 
  int height = 700;
  int width = 700;
  double xmin = 0.0;
  double xmax = 2.0;
  double ymin = 0.0;
  double ymax = 2.0;
  double zmin;
  double zmax;
  double xw, yh, xoff, yoff;
  int nx = 9;
  int ny = 9;
  double hx, hy;
  double ya = 0.5; //  orthographic typ 0.2 to 0.7, must fix text 
  double xs = (0.666*(width-100));
  double ys = (0.666*(height-100));
  int xp, yp; // for passing back ortho point

  public test_bihar2d()
  {
    System.out.println("test_bihar2d.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, ux and uy, the exterior normal derivatives");
    System.out.println("   of u must be specified at the boundary.");
    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);

    setTitle("test_bihar2d");
    setSize(width,height);
    setBackground(Color.white);
    setForeground(Color.black);
    addWindowListener(new WindowAdapter()
    {
      public void windowClosing(WindowEvent e)
      {
        System.exit(0);
      }
    });
    setVisible(true);
    this.addMouseListener (new mousePressHandler());

    do_program();
  } // test_bihar2d constructor

  void do_program()
  {
    double x, y, val, err, maxerr, avgerr, rmserr;
    double xp, yp , sumerr, sumsqerr, deriv, eval_deriv;
    int i_max, j_max, cnt;

    hx = (xmax-xmin)/(double)(nx-1);
    hy = (ymax-ymin)/(double)(ny-1);
    System.out.println("nx="+nx+", ny="+ny+", hx="+hx+", hy="+hy);
    // for plotting
    xw = (xmax-xmin)*1.2;
    yh = (ymax-ymin)*1.2;
    xoff = 0.083333*xw;
    yoff = 0.083333*yh;
    // for derivative check
    double xg[] = new double[nx];
    double yg[] = new double[ny];
    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 xxyy
    double ug[][] = new double[nx][ny];

    zmin = 100.0;
    zmax = -100.0;
    // check that analytic solution numerically satisfies PDE
    maxerr = 0.0;
    avgerr = 0.0;
    i_max = 0;
    j_max = 0;
    for(int i=0; i<nx; i++)
    {
      x = xmin + (double)i*hx;
      xg[i] = x;
      for(int j=0; j<ny; j++)
      {
	y = ymin + (double)j*hy;
        yg[j] = y;
        err = uxxxx(x,y) + 2.0*uxxyy(x,y) + uyyyy(x,y) + 
	      2.0*uxx(x,y) + 2.0*uyy(x,y)  + 2.0*u(x,y) - f(x,y);
        ug[i][j] = ub(x,y);
        err = Math.abs(err); // PDE check against f(x,y)
        avgerr += err;
	if(err>maxerr) 
	{
	  maxerr = err;
	  i_max = i;
	  j_max = j;
	}
        zmax = Math.max(zmax,u(x,y));
        zmin = Math.min(zmin,u(x,y));
      }
    }
    x = xmin+i_max*hx;
    y = ymin+j_max*hy;
    System.out.println("maxerr PDE-f(x,y) is "+maxerr+", at i_max="+
                       i_max+", j_max="+j_max);
    System.out.println("x="+x+", y="+y+", f(x,y)="+f(x,y));
    System.out.println("zmin="+zmin+", zmax="+zmax+", avgerr="+(avgerr/(nx*ny)));

    // check numerical accuracy of computation of derivatives
    sumerr = 0.0;
    sumsqerr = 0.0;
    cnt = 0;
    for(int i=1; i<nx-1; i++)
    {
      // compute derivative coefficients in x direction
      rderiv(2,nx,i,hx,cxx[i]);
      // for(int k=0; k<nx; k++)
      //   System.out.println("cxx point i="+i+", term k="+k+", coef="+cxx[i][k]); 
      rderiv(4,nx,i,hx,cxxxx[i]);
    }
    for(int ii=1; ii<ny-1; ii++)
    {
      // compute derivative coefficients in y direction
      rderiv(2,ny,ii,hy,cyy[ii]);
      rderiv(4,ny,ii,hy,cyyyy[ii]);
    }

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uxx
    {
      for(int j=1; j<ny-1; j++)
      {
	eval_deriv = 0.0; // numeric
        deriv = uxx(xg[i],yg[j]); // analytic
        for(int k=0; k<nx; k++)
        {
	  eval_deriv += cxx[i][k]*ub(xg[k],yg[j]);
	}
        err = deriv - eval_deriv;
        maxerr = Math.max(maxerr,Math.abs(err));
      }
    }
    System.out.println("check cxx  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uxxxx
    {
      for(int j=1; j<ny-1; j++)
      {
	eval_deriv = 0.0; // numeric
        deriv = uxxxx(xg[i],yg[j]); // analytic
        for(int k=0; k<nx; k++)
        {
	  eval_deriv += cxxxx[i][k]*ub(xg[k],yg[j]);
	}
        err = deriv - eval_deriv;
        maxerr = Math.max(maxerr,Math.abs(err));
      }
    }
    System.out.println("check cxxxx  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uyy
    {
      for(int j=1; j<ny-1; j++)
      {
	eval_deriv = 0.0; // numeric
        deriv = uyy(xg[i],yg[j]); // analytic
        for(int k=0; k<ny; k++)
        {
	  eval_deriv += cyy[j][k]*ub(xg[i],yg[k]);
	}
        err = deriv - eval_deriv;
        maxerr = Math.max(maxerr,Math.abs(err));
      }
    }
    System.out.println("check cyy  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uyyyy
    {
      for(int j=1; j<ny-1; j++)
      {
	eval_deriv = 0.0; // numeric
        deriv = uyyyy(xg[i],yg[j]); // analytic
        for(int k=0; k<ny; k++)
        {
	  eval_deriv += cyyyy[j][k]*ub(xg[i],yg[k]);
	}
        err = deriv - eval_deriv;
        maxerr = Math.max(maxerr,Math.abs(err));
      }
    }
    System.out.println("check cyyyy  maxerr="+maxerr);

    maxerr = 0.0;
    for(int i=1; i<nx-1; i++) // check uxxyy
    {
      for(int j=1; j<ny-1; j++)
      {
        deriv = uxxyy(xg[i],yg[j]); // analytic
        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]*ub(xg[k],yg[jj]);
	  }
          if(Math.abs(uxx(xg[i],yg[jj])-tcxxyy[jj])>0.001)
	    System.out.println("tcxxy["+jj+"]="+tcxxyy[jj]+
			       " != uxx="+uxx(xg[i],yg[jj]));
	}
        eval_deriv = 0.0;
        for(int k=0; k<ny; k++)  // uyy part of uxxyy
        {
	  eval_deriv += cyy[j][k]*tcxxyy[k];
	}
        err = deriv - eval_deriv;
        maxerr = Math.max(maxerr,Math.abs(err));
      }
    }
    System.out.println("check cxxyy  maxerr="+maxerr);

    for(int i=1; i<nx-1; i++) // full equation
    {
      x = xmin + (double)i*hx;
      for(int j=1; j<ny-1; j++)
      {
	y = ymin + (double)j*hy;
	eval_deriv = 2.0*ug[i][j]; // PDE u term
        
        for(int k=0; k<nx; k++)
        {
	  eval_deriv += 2.0*cxx[i][k]*ug[k][j]; // PDE  uxx term
	  eval_deriv += cxxxx[i][k]*ug[k][j];   // PDE  uxxxx term
	} 
        for(int k=0; k<ny; k++)
        {
	  eval_deriv += 2.0*cyy[j][k]*ug[i][k];  // PDE uyy term
	  eval_deriv += cyyyy[j][k]*ug[i][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][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(xg[i],yg[j]);
        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 all terms together in this PDE");
    System.out.println("maxerr="+maxerr+", rmserr="+rmserr+", avgerr="+avgerr);

    // plot
    repaint();

  } // end do_program

  class mousePressHandler extends MouseAdapter
  {
    public void mousePressed (MouseEvent e)
    {
      xm = e.getX();
      ym = e.getY();
      bm = e.getButton();
      // System.out.println("down xm="+xm+", ym="+ym+", bm="+bm); // debug

      if(in_rect(width-70, height-20, 30, 20))
      {
	// do something, run
	requestFocus();
        repaint();
      }
      else if(in_rect(width-105, height-20, 30, 20))
      {
	// do something, next
	requestFocus();
	nx = nx+2;
	ny = ny+2;
        do_program();
        repaint();
      }
      else if(in_rect(width-145, height-20, 35, 20))
      {
	// do something, back
	requestFocus();
        if(nx>5)
	{
	  nx = nx-2;
	  ny = ny-2;
	}
        do_program();
        repaint();
      }
      if(in_rect(width-35, height-20, 30, 20)) System.exit(0);
    } // end mousePressed 
  } // end mousePressHandler

  boolean in_rect(int x, int y, int w, int h)
  {
    return xm>x && xm<x+w && ym>y && ym<y+h;
  } // end in_rect

  boolean not_in(int i, int a[], int n)
  {
    boolean found = false;
    int j;

    for(j=0; j<n; j++)
    {
      if(i==a[j]) { found=true; break; }
    }
    return !found;
  } // end not_in

  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;
  }

  double u(double x, double y)
  {
    return ub(x,y);
  }

  double ux(double x, double y)
  {
    return 0.04*x*x*x - 0.08*x*y*y + 0.3 *x*x*y + 0.2*y*y*y - y;
  }

  double uxx(double x, double y)
  {
    return 0.12*x*x -0.08*y*y + 0.6*x*y;
  }

  double uxxx(double x, double y)
  {
    return 0.24*x + 0.6*y;
  }

  double uxxxx(double x, double y)
  {
    return 0.24;
  }

  double uy(double x, double y)
  {
    return 0.08*y*y*y - 0.8*x*x*y + 0.1*x*x*x + 0.6*x*y*y -x;
  }

  double uyy(double x, double y)
  {
    return 0.24*y*y -0.08*x*x + 1.2*x*y;
  }

  double uyyy(double x, double y)
  {
    return 0.48*y + 1.2*x;
  }

  double uyyyy(double x, double y)
  {
    return 0.48;
  }

  double uxxyy(double x, double y)
  {
    return -0.16;
  }
 
  public void paint(Graphics g)
  {
    double x, y, z; // z here is "u"
    double z1, z2, hz;
    int xp1, yp1, xp2, yp2; // line segment
    hx = (xmax-xmin)/nx; // for plotting with x+hx
    hy = (ymax-ymin)/ny;

    // buttons
    g.setColor(Color.red);
    g.fillRect(width-35,height-20,30,15);
    g.setColor(Color.black);
    g.drawString("quit", width-30, height-10);
    g.setColor(Color.green);
    g.fillRect(width-70,height-20,30,15);
    g.setColor(Color.black);
    g.drawString("run", width-65, height-10);
    g.setColor(Color.yellow);
    g.fillRect(width-105,height-20,30,15);
    g.setColor(Color.black);
    g.drawString("fine", width-100, height-10);
    g.setColor(Color.orange);
    g.fillRect(width-145,height-20,35,15);
    g.setColor(Color.black);
    g.drawString("big ", width-140, height-10);

    // axis labels
    g.drawString("xmin", 25, height-15);
    g.drawString("xmax", (int)(0.66*width-70), height-15);
    g.drawString("ymin", (int)(0.66*width-30), height-25);
    g.drawString("ymax", (int)(width-80), (int)(0.66*height+25));
    g.drawString("umin", (int)(width-75), (int)(0.66*height+5));
    g.drawString("umax", (int)(width-75), (int)(85));

    // height color scale
    int hgt = 210;
    int ncol = 10;
    for(int i=0; i<ncol; i++)
    {
      set_color(i, g);
      g.fillRect(20,hgt-20*i,50,15);
    }
    g.drawString("u value", 80, 40);

    // draw box, bottom
    set_color(0, g);
    ortho(xmin, ymin, zmin); xp1=xp; yp1=yp;
    ortho(xmax, ymin, zmin); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmax, ymax, zmin); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    ortho(xmin, ymax, zmin); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmin, ymin, zmin); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    set_color(ncol-1, g); // top
    ortho(xmin, ymin, zmax); xp1=xp; yp1=yp;
    ortho(xmax, ymin, zmax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmax, ymax, zmax); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    ortho(xmin, ymax, zmax); xp2=xp; yp2=yp;
    g.drawLine(xp1, yp1, xp2, yp2);
    ortho(xmin, ymin, zmax); xp1=xp; yp1=yp;
    g.drawLine(xp2, yp2, xp1, yp1);
    // draw box edges with z value colors
    z1 = zmin;
    hz = (zmax-zmin)/(ncol);
    for(int i=0; i<ncol; i++)
    {
      z2 = z1+hz;
      set_color(i,g);
      ortho(xmin, ymin, z1); xp1=xp; yp1=yp;
      ortho(xmin, ymin, z2); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      ortho(xmin, ymax, z1); xp1=xp; yp1=yp;
      ortho(xmin, ymax, z2); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      ortho(xmax, ymax, z1); xp1=xp; yp1=yp;
      ortho(xmax, ymax, z2); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      ortho(xmax, ymin, z1); xp1=xp; yp1=yp;
      ortho(xmax, ymin, z2); xp2=xp; yp2=yp;
      g.drawLine(xp1, yp1, xp2, yp2);
      z1 = z2;
    }

    // plot function or data with z value colors
    for(int i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      for(int j=0; j<ny; j++)
      {
	y = ymin+(double)j*hy;
        z1 = ub(x, y);
        ortho(x, y, z1); xp1=xp; yp1=yp;
        z2 = ub(x+hx, y);
        ortho(x+hx, y, z2); xp2=xp; yp2=yp;
        set_color((int)((ncol*(((z1+z2)/2.0)-zmin)/(zmax-zmin))+0.5), g);
        g.drawLine(xp1, yp1, xp2, yp2);
        z2 = ub(x, y+hy);
        ortho(x, y+hy, z2); xp2=xp; yp2=yp;
        set_color((int)((ncol*(((z1+z2)/2.0)-zmin)/(zmax-zmin))+0.5), g);
        g.drawLine(xp1, yp1, xp2, yp2);
      }
    }
    for(int i=0; i<nx; i++)
    {
      x = xmin+(double)i*hx;
      y = ymax;
      z1 = ub(x, y);
      ortho(x, y, z1); xp1=xp; yp1=yp;
      z2 = ub(x+hx, y);
      ortho(x+hx, y, z2); xp2=xp; yp2=yp;
      set_color((int)((ncol*(((z1+z2)/2.0)-zmin)/(zmax-zmin))+0.5), g);
      g.drawLine(xp1, yp1, xp2, yp2);
    }
    for(int j=0; j<ny; j++)
    {
      y = ymin+(double)j*hy;
      x = xmax;
      z1 = ub(x, y);
      ortho(x, y, z1); xp1=xp; yp1=yp;
      z2 = ub(x, y+hy);
      ortho(x, y+hy, z2); xp2=xp; yp2=yp;
      set_color((int)((ncol*(((z1+z2)/2.0)-zmin)/(zmax-zmin))+0.5), g);
      g.drawLine(xp1, yp1, xp2, yp2);
    }
  } // end paint

  // orthographic projection, sets xp,yp
  void ortho(double x, double y, double z)
  {
    double x1, y1, z1, x2, y2, x3, y3;
    x1 = (x-xmin)/(xmax-xmin);   // zero base, normalized
    y1 = (y-ymin)/(ymax-ymin); 
    z1 = (z-zmin)/(zmax-zmin);  // z-zmin;
    x2 = x1+y1*ya; // ortho z is up, y is back
    y2 = z1+y1*ya;
    x3 = x2*xs;    // scale 
    y3 = y2*ys;
    x3 = x3+25;        // away from boarder
    y3 = height-y3-25; // away from boarder
    xp = (int)x3;
    yp = (int)y3;
    //System.out.println("ortho xp="+xp+", yp="+yp);
    //System.out.println("x1="+x1+", y1="+y1+", z1="+z1);
    //System.out.println(" ");
  } // end ortho

  void set_color(int i, Graphics g)
  {
    if(i<0) {g.setColor(Color.black); return;}
    if(i>9) {g.setColor(Color.gray); return; }
    switch(i)
    {
      case 0: g.setColor(Color.black); break;
      case 1: g.setColor(Color.pink); break; // should be brown
      case 2: g.setColor(Color.red); break;
      case 3: g.setColor(Color.orange); break;
      case 4: g.setColor(Color.yellow); break;
      case 5: g.setColor(Color.green); break;
      case 6: g.setColor(Color.blue); break;
      case 7: g.setColor(Color.magenta); break;
      case 8: g.setColor(Color.cyan); break;
      case 9: g.setColor(Color.gray); break;
    }
  } // end set_color

  void rderiv(int order, int npts, int point, double h, double c[])
  {
    nderiv D = new nderiv();
    int a[] = new int[npts]; // basic discretization from nderiv
    int b[] = new int[1];    // out parameter, denominator

    D.deriv(order,npts,point,a,b);
    for(int i=0; i<npts; i++)
    {
      c[i] = (double)a[i]/((double)b[0]);
      for(int j=0; j<order; j++) c[i] = c[i]/h;
    }
  } // end rderiv
 
  public static void main(String[] args)
  {
    new test_bihar2d();
  }
} // end class  test_bihar2d.java