// w0015pde.java  airfoil in wind tunnel, non uniform grid
//                 data read from w0015a22.grid produced by w0015grid.java
//                 edit  ndel=10 e.g.  order=5 e.g
//                 airfoil has a chord of 1 meter,

import java.awt.*;
import java.awt.event.*;
import java.io.*;

// needs invert.java  nuderiv.java  simeq_newton5.java  read_grid.java
// and a data file  w0015a22.grid

public class w0015pde extends Frame
{
  double xmin = -1.0;
  double xmax = 2.0;
  double xoff = 340.0;
  double yoff = 500.0;
  double ymin = -1.25;
  double ymax = 1.15;
  double scale = 320.0; // for plotting
  double xminw; // from data read
  double xmaxw;
  double yminw;
  double ymaxw;
  double xadj = -0.5; // approx center for output 
  double grid = 0.01; // grid spacing, may be computed 
  int width = 1000;
  int height = 1000; // upper part for title 

  String fname; // .grid file has the folloairfoil
  int ixbound[] = new int[100]; // output boundary
  int jybound[] = new int[100];
  double xg[] = new double[100]; // output grid
  double yg[] = new double[100];
  int nxg;
  int nyg;
  int nbound;

  int ndel = 4; // grids outside airfoil
  double angle_deg = 22.0; // may be read in 
  String angle_ch = "22";
  int nDOF; // number of degrees of freedom, u,v,p independent variable each
  int neqn; // number of equations 3*nDOF = number of variables
  int nlin; // number of nonlinear terms u_i*v_j for some i,j, some uii, vjj
  int ntot; // independent plus product variables neqn+ntot

  int order = 5; // 5 fifth order nonuniform solution, some code depends
  double cx[];   // for nuderiv
  double cxx[];  // computed by discrete
  double cy[];
  double cyy[];
  double xx[];
  double yy[];
  int iu[];      // x index is u position
  int jv[];      // y index+nDOF is v position
  int eqn = 0;   // equation being built
  int var1[];    // variable number for non linear equations
  int var2[];    // variable number for non linear equations
  int var3[];    // variable number for non linear equations
  int vari1[];   // variable number for non linear equations
  int vari2[];   // variable number for non linear equations
  double Y[];    // RHS right hand side of equations
  double X[];    // initial guess and solution of non linear equations
  double A[][];  // A Y = X  system of non linear equations
  int nl;        // available non linear term in var1, var2
                 // v is offset, nDOF+j
  double u0 = 88.0;     // u velocity, uniform at input, output, top, bottom
  double v0 = 0.0;      // v velocity is zero at these boundaries
  double rho = 1.223;   // density
  double mu = 182.0E-6; // viscosity
  double p0 = 0.100;    // pressure at input and output, Pa about 14.5 psi
  simeq_newton5 SN = new simeq_newton5(); // may change parameters 

  w0015pde(String gname, int tndel, int torder)
  {
    System.out.println("w0015pde.java running");
    fname = gname;
    ndel = tndel;
    order = torder;
    System.out.println("using "+fname+", ndel="+ndel+", order="+order);
    cx = new double[order];   // for nuderiv
    cxx = new double[order];  // computed by discrete
    cy = new double[order];
    cyy = new double[order];
    xx = new double[order];
    yy = new double[order];
    iu = new int[order];      // x index is u position
    jv = new int[order];      // y index+nDOF is v position
    build_vert();  
    build_equations();

    setTitle(fname+" "); // ndel, order  to String
    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());

  } // end constructor w0015pde

  void build_vert()
  {
    double x, y;

    // read in a airfoil
    System.out.println("about to read airfoil grid from file "+fname);
    read_grid G = new read_grid(fname);
    nxg = G.read_xgrid(xg);
    nyg = G.read_ygrid(yg);
    nbound = G.read_bound(ixbound, jybound);
    for(int k=0; k<nbound; k++) System.out.println("ixb["+k+"]="+ixbound[k]+
                                                   ", jyb="+jybound[k]);

    // offset and build surrounding wind tunnel grid
    if(nxg>0)
    {
      xminw = xg[0];
      xmaxw = xg[nxg-1];
      yminw = yg[0];
      ymaxw = yg[nyg-1];
    }
    else
    {
      xminw = 0.0; // null grid
      xmaxw = 0.0;
      yminw = 0.0;
      ymaxw = 0.0;
    }
    System.out.println("xminw="+xminw+", xmaxw="+xmaxw+
                       ", yminw="+yminw+", ymaxw="+ymaxw);
    // xmin, xmax, ymin, ymax set manually for now
    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax);
    double xmind = -(xmin-xminw)/(double)(ndel);
    double xmaxd = (xmax-xmaxw)/(double)(ndel);
    double ymind = -(ymin-yminw)/(double)(ndel);
    double ymaxd = (ymax-ymaxw)/(double)(ndel);
    System.out.println("xmind="+xmind+", xmaxd="+xmaxd+
                       ", ymind="+ymind+", ymaxd="+ymaxd);
    // move airfoil and insert additional grid
    for(int i=nxg-1; i>=0; i--) xg[i+ndel] = xg[i]; 
    for(int i=nyg-1; i>=0; i--) yg[i+ndel] = yg[i]; 
    for(int i=0; i<nbound; i++) ixbound[i] = ixbound[i]+ndel;
    for(int i=0; i<nbound; i++) jybound[i] = jybound[i]+ndel;
    x = xmin;
    for(int i=0; i<ndel; i++) {xg[i]=x; x+=xmind;};
    x = xmaxw+xmaxd;
    for(int i=ndel+nxg; i<ndel+ndel+nxg; i++) {xg[i]=x; x+=xmaxd;};
    nxg = nxg+ndel+ndel;
    for(int i=0; i<nxg; i++) System.out.println("xg["+i+"]="+xg[i]);

    y = ymin;
    for(int i=0; i<ndel; i++) {yg[i]=y; y+=ymind;};
    y = ymaxw+ymaxd;
    for(int i=ndel+nyg; i<ndel+ndel+nyg; i++) {yg[i]=y; y+=ymaxd;};
    nyg = nyg+ndel+ndel;
    for(int j=0; j<nyg; j++) System.out.println("yg["+j+"]="+yg[j]);
  } // end build_vert
    
  int bound(int i, int j) // 0 for DOF, 1 for external, 2 for airfoil, 3 inside
  {
    int jval;

    if(i==0 || i==nxg-1) return 1;
    if(j==0 || j==nyg-1) return 1;
    if(i<ndel || j<ndel) return 0;
    if(i>nxg-ndel || j>nyg-ndel) return 0;
    for(int k=0; k<nbound; k++)
      if(i==ixbound[k] && j==jybound[k]) return 2;
    jval=0;
    for(int k=0; k<nbound; k++)
    {
      if(i==ixbound[k])
      {
	 if(jval==0) jval = jybound[k];
         else
	 {
	   // may not be any
           if((j<jval&&j>jybound[k]) || (j<jybound[k]&&j>jval)) return 3;
           jval = 0;
	 }
      }
    }
    return 0;
  } // end bound

  //  wind tunnel PDE, Navier Stokes static, incompressible, 2D
  //            area defined by xg[0]..xg[nx-1]  yg[0]..yg[ny-1]
  //            at every interior  i,j  xg[i],yg[j]  x,y :
  //            eq1 rho(u(x,y)*ux(x,y) + v(x,y)*uy(x,y)) + px(x,y) - 
  //                mu*(uxx(x,y) +uyy(x,y)) = 0
  //            eq2 rho(u(x,y)*vx(x,y) + v(x,y)*vy(x,y)) + py(x,y) -
  //                mu*(vxx(x,y) +vyy(x,y)) = 0
  //            eq3 ux(x,y) + vy(x,y) = 0   continuity equation
  //
  // u(x,y) is velocity in meters/second in X direction at  x,y
  // ux(x,y) is first derivative of u(x,y) wrt x
  // uy(x,y) is first derivative of u(x,y) wrt y
  // uxx(x,y) is second derivative of u(x,y) wrt x
  // uyy(x,y) is second derivative of u(x,y) wrt y
  //
  // v(x,y) is velocity in meters/second in Y direction at  x,y
  // vx(x,y) is first derivative of v(x,y) wrt x
  // vy(x,y) is first derivative of v(x,y) wrt y
  // vxx(x,y) is second derivative of v(x,y) wrt x
  // vyy(x,y) is second derivative of v(x,y) wrt y
  //
  // p(x,y) is relative pressure at  x,y in kilogram/(meter second^2)
  // px(x,y) is first derivative of p(x,y) wrt x
  // py(x,y) is first derivative of p(x,y) wrt y
  //
  // rho is air density   1.223 kilogram/meter^3  at 15.7 degrees Celsius
  // mu  is air viscosity 182*10^-6 kilogram/(meter second)
  // u0 input and output X velocity 88 meters per second,
  // v0 all boundary Y velocity 0.0 meters per second
  // p0 external pressure 100kPa about 14.5 psi  0.1 Kg/meter^2
  //
  // force of gravity is ignored
  // change in density rho due to pressure, ignored in first cut.
  // later, use gas law and iterate 

  void build_equations()
  {
    double x, y;
    int ixpt[] = new int[nxg*nyg]; // max size, actual nDOF
    int jypt[] = new int[nxg*nyg];

    nDOF = 0;
    for(int i=0; i<nxg; i++)
    {
      for(int j=0; j<nyg; j++)
      {
        if(i==0 && j== 0)
           System.out.println("in build, x="+((int)(scale*xg[i]+xoff))+
                              ", y="+(height-(int)(scale*yg[j]+yoff)));
        if(bound(i,j)!=0) continue; // not a DOF
        // if(bound(i,j)==1) outer boundary
        // if(bound(i,j)==2) airfoil edge
        // if(bound(i,j)==3) hole
        x = xg[i];
        y = yg[j];
        ixpt[nDOF] = i;
        jypt[nDOF] = j;

        nDOF++;
      }
    }
    System.out.println("DOF = "+nDOF);
    // for Navier Stokes in two dimensions, u, v, p at each DOF
    // for Navier Stokes three equations for each DOF
    neqn = 3*nDOF; // also number of linear unknowns
    nlin = nDOF*nDOF; // actual maximum less than order*order*nDOF;
    ntot = neqn + nlin; // maximum columnc with order u^2 v^2 u*v terms
    System.out.println("neqn="+neqn+", ntot="+ntot+", A is "+
                       (neqn*ntot*8)+" bytes");
    var1 = new int[ntot];
    var2 = new int[ntot];
    var3 = new int[ntot];
    vari1 = new int[ntot];
    vari2 = new int[ntot];
    Y = new double[neqn];
    X = new double[ntot];
    A = new double[neqn][ntot];
    nl = neqn; // available non linear term in var1, var2
               // v is offset, nDOF+j

    // set up variable and product definition in var1 and var2
    // nDOF u, nDOF v, nDOF p,
    // startine with [neqn] u_0*v_0 u_0*v_1 ... u_0*v_nDOF as required 
    for(int k=0; k<nDOF; k++)
    {
      var1[k] = k; // u
      var2[k] = -1;
      var3[k] = -1;
      vari1[k] = -1;
      vari2[k] = -1;
      var1[k+nDOF] = k+nDOF; // v
      var2[k+nDOF] = -1;
      var3[k+nDOF] = -1;
      vari1[k+nDOF] = -1;
      vari2[k+nDOF] = -1;
      var1[k+2*nDOF] = k+2*nDOF; // p
      var2[k+2*nDOF] = -1;
      var3[k+2*nDOF] = -1;
      vari1[k+2*nDOF] = -1;
      vari2[k+2*nDOF] = -1;
    }
    nl = neqn;
    for(int i=0; i<nDOF; i++)
    {
      for(int j=0; j<nDOF; j++)
      {
        var1[nl] = i; // must be set for nonlinear term
        var2[nl] = j; // must be set for nonlinear term 
        var3[nl] = -1; // most A[i][j] zero, ecomomize later
	nl++;
      }
    }

    for(int i=1; i<nxg-1; i++) // generate equations, 3, for each DOF
    {                          // skip boundaries and internal
      for(int j=1; j<nyg-1; j++)
      {
        if(bound(i,j)!=0) continue; // not a DOF
        // if(bound(i,j)==1) outer boundary
        // if(bound(i,j)==2) airfoil edge
        // if(bound(i,j)==3) hole
        discrete(i,j); // fill cx, cxx, cy, cyy, xx, yy, iu, jv
        if((i==1 && j==1) || (i==2 && j==1) || (i==1 && j==nyg-2) ||
           (i==nxg-2 && j==1) || (i==nxg-2 && j==nyg-2))
	{
	  System.out.println(" x iu     xx      cx       cxx");
	  System.out.println(" y jv     yy      cy       cyy");
          for(int k=0; k<order; k++)
	  {
	    System.out.println("x "+iu[k]+" "+xx[k]+" "+cx[k]+" "+cxx[k]);
	    System.out.println("y "+jv[k]+" "+yy[k]+" "+cy[k]+" "+cyy[k]);
	  }
	} // end if for debug print

        fill_matrix(i, j); // eq1, eq2, eq3 A[][] and Y is RHS
      }
    }

    // initialize X  all u=u0,  all v=v0  all p=p0
    //               use var1 and var2 to initialize products
    for(int k=0; k<nDOF; k++)
    {
      X[k] = u0;
      X[k+nDOF] = v0;
      X[k+2*nDOF] = p0;
    }
    for(int k=neqn; k<ntot; k++)
    {
      if(var1[k]>=0)
      {
	X[k] = X[var1[k]]*X[var2[k]];
      }
    }

    // solve nonlinear equations, use results u, v, p
    SN.simeq(neqn, ntot-neqn, A, Y, var1, var2, var3, vari1, vari2, X);
    //print results from X
    System.out.println(" ix  jy    u           v          p");
    for(int i=0; i<nDOF; i++)
    {
      System.out.println(ixpt[i]+"  "+jypt[i]+"  "+X[i]+" "+
                         X[i+nDOF]+" "+X[i+2*nDOF]);
    }
  } // end build_equations

  void discrete(int ix, int jy) // xg,yg indices
  {
    int ixu = ix*nyg+jy;      // linear pointer in X solution
    int jyv = ix*nyg+jy+nDOF; // linear pointer in X solution
    // try for ix,jy in center  order points (5 for now)
    if(bound(ix-1,jy)==0) // x before OK
    {
      if(bound(ix+1,jy)==0) // x after OK
      {
	for(int i=0; i<order; i++)
	{
	  iu[i] = ixu+i-2;
          xx[i] = xg[i+ix-2];
	}
        new nuderiv(1, order, 2, xx, cx);
        new nuderiv(2, order, 2, xx, cxx);
      }
      else // x after is boundary, move back 1
      {
	if(bound(ix-2,jy)!=0)
          System.out.println("ERROR discrete x at ix="+ix+", jy="+jy);
	for(int i=0; i<order; i++)
	{
	  iu[i] = ixu+i-3;
          xx[i] = xg[i+ix-3];
	}
        new nuderiv(1, order, 1, xx, cx);
        new nuderiv(2, order, 1, xx, cxx);
      }
    }
    else // x before is boundary, move forward one
    {
      if(bound(ix+2,jy)!=0)
        System.out.println("ERROR discrete x at ix="+ix+", jy="+jy);
      for(int i=0; i<order; i++)
      {
        iu[i] = ixu+i-1;
        xx[i] = xg[i+ix-1];
      }
      new nuderiv(1, order, 3, xx, cx);
      new nuderiv(2, order, 3, xx, cxx);
    }
    if(bound(ix,jy-1)==0) // y before OK
    {
      if(bound(ix,jy+1)==0) // y after OK
      {
	for(int j=0; j<order; j++)
	{
	  jv[j] = jyv+j-2;
          yy[j] = yg[j+jy-2];
	}
        new nuderiv(1, order, 2, yy, cy);
        new nuderiv(2, order, 2, yy, cyy);
      }
      else // y after is boundary, move back 1
      {
	if(bound(ix,jy-2)!=0)
          System.out.println("ERROR discrete at y ix="+ix+", jy="+jy);
	for(int j=0; j<order; j++)
	{
	  jv[j] = jyv+j-3;
          yy[j] = yg[j+jy-3];
	}
        new nuderiv(1, order, 1, yy, cy);
        new nuderiv(2, order, 1, yy, cyy);
      }
    }
    else // y before is boundary, move forward one
    {
      if(bound(ix,jy+2)!=0)
        System.out.println("ERROR discrete y at ix="+ix+", jy="+jy);
      for(int j=0; j<order; j++)
      {
        jv[j] = jyv+j-1;
        yy[j] = yg[j+jy-1];
      }
      new nuderiv(1, order, 3, yy, cy);
      new nuderiv(2, order, 3, yy, cyy);
    }
  } // end discrete

  void fill_matrix(int ix, int jy)
  {
    int ixu = ix*nyg+jy;  // linear subscript of u in X solution
    int jyv = ixu+nDOF;   // linear subscript of v in X solution
    int kp  = ixu+2*nDOF; // linear subscript of p in X solution

    // eq1:  +rho * u(x,y)*ux(x,y)  cx  check var1, var2, add if new
    //       +rho * v(x,y)*uy(x,y)  cy  "
    //       +px(x,y)               cx
    for(int k=0; k<order; k++) A[eqn][iu[k]] +=cx[k]; // ? kp
    //       -mu  * uxx(x,y)        cxx
    if(iu[0]==0) // boundary
    {
      Y[eqn] += u0*mu*cxx[0];
      for(int k=1; k<order; k++) A[eqn][iu[k]] -= mu*cxx[k];
    }
    else if(iu[order-1]==nxg-1)
    {
      Y[eqn] += u0*mu*cxx[order-1];
      for(int k=0; k<order-1; k++) A[eqn][iu[k]] -= mu*cxx[k];
    }
    for(int k=0; k<order; k++) A[eqn][iu[k]] -= mu*cxx[k];
    //       -mu  * uyy(x,y))       cyy
    if(iu[0]==0) // boundary
    {
      Y[eqn] += u0*mu*cyy[0];
      for(int k=1; k<order; k++) A[eqn][iu[k]] -= mu*cyy[k];
    }
    else if(iu[order-1]==nxg-1)
    {
      Y[eqn] += u0*mu*cyy[order-1];
      for(int k=0; k<order-1; k++) A[eqn][iu[k]] -= mu*cyy[k];
    }
    for(int k=0; k<order; k++) A[eqn][iu[k]] -= mu*cyy[k];
    //       RHS = 0 unless boundary
    eqn++;
    // eq2:  +rho * u(x,y)*vx(x,y)  cx  "
    //       +rho * v(x,y)*vy(x,y)  cy  "
    //       +py(x,y)               cy
    for(int k=0; k<order; k++) A[eqn][jv[k]] +=cy[k]; // ? kp
    //       -mu  * vxx(x,y)        cxx
    if(jv[0]==0) // boundary
    {
      Y[eqn] += mu*v0*cxx[0];
      for(int k=1; k<order; k++) A[eqn][jv[k]] -= mu*cxx[k];
    }
    else if(jv[order-1]==nyg-1)
    {
      Y[eqn] += mu*v0*cxx[order-1];
      for(int k=0; k<order-1; k++) A[eqn][jv[k]] -= mu*cxx[k];
    }
    for(int k=0; k<order; k++) A[eqn][jv[k]] -= mu*cxx[k];
    //       -mu  * vyy(x,y)        cyy
    for(int k=0; k<order; k++) A[eqn][jv[k]] -=(mu*cyy[k]);
    //        RHS = 0 unless boundary
    eqn++;
    // eq3:   +ux(x,y)               cx
    if(iu[0]==0) // boundary
    {
      Y[eqn] -= u0*cx[0];
      for(int k=1; k<order; k++) A[eqn][iu[k]] +=cx[k];
    }
    else if(iu[order-1]==nxg-1)
    {
      Y[eqn] -= u0*cx[order-1];
      for(int k=0; k<order-1; k++) A[eqn][iu[k]] +=cx[k];
    }
    for(int k=0; k<order; k++) A[eqn][iu[k]] +=cx[k];
    //        +vy(x,y)               cy
    if(jv[0]==0) // boundary
    {
      Y[eqn] -= v0*cy[0];
      for(int k=1; k<order; k++) A[eqn][jv[k]] +=cy[k];
    }
    else if(jv[order-1]==nyg-1)
    {
      Y[eqn] -= v0*cy[order-1];
      for(int k=0; k<order-1; k++) A[eqn][jv[k]] +=cy[k];
    }
    for(int k=0; k<order; k++) A[eqn][jv[k]] +=cy[k];
    //        RHS = 0 unless boundary
    eqn++;

  } // end fill_matrix

  class mousePressHandler extends MouseAdapter
  {
    public void mousePressed (MouseEvent e)
    {
      int x, y, b;
      x = e.getX();
      y = e.getY();
      b = e.getButton();
      if(b==3)
      {
      }
      requestFocus();
      System.out.println("at x="+x+"   y="+y+"   b="+b); // debug print
      repaint();
    } // end mousePressed
  } // end class mousePressedHandler

  public void paint(Graphics g)
  {
    int ix,jy;
    g.setColor(Color.black);
    g.drawString(fname, 20, 50);
    g.setColor(Color.red);
    for(int i=0; i<nxg; i++)
    {
      for(int j=0; j<nyg; j++)
      {
        if(i==0 && j== 0)
           System.out.println("in paint, x="+((int)(scale*xg[i]+xoff))+
                              ", y="+(height-(int)(scale*yg[j]+yoff)));
        if(bound(i,j)==0) g.setColor(Color.red);
        if(bound(i,j)==1) g.setColor(Color.black);
        if(bound(i,j)==2) g.setColor(Color.blue);
        if(bound(i,j)==3) g.setColor(Color.yellow);
 
        g.drawOval((int)(scale*xg[i]+xoff), 
	  	   height-(int)(scale*yg[j]+yoff),
                   1, 1);
      } // end j
    } // end i
  } // end paint

  public static void main(String args[])
  {
    // read in name of .grid file
    String gname = "w0015a22.grid";
    if(args.length>0)
        gname = args[0];
    int tndel = 4; // minimum number of points from airfoil to boundary 
    if(args.length>1)
	tndel = Integer.parseInt(args[1]);
    int torder = 5; // discretion order
    if(args.length>2)
	torder = Integer.parseInt(args[2]);
    if(torder>tndel+1) torder = tndel+1;

    new w0015pde(gname, tndel, torder);
    System.out.println("w0015pde.java ends");
  } // end main
} // end class w0015pde