// cylinder_flow.java  cylinder in wind tunnel, uniform grid
//                 hard coded circle, two D view of cylinder

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

public class cylinder_flow extends Frame
{
  int order = 7; // 7th order approximation of derivatives
  double cx[]  = new double[order];   // for rderiv
  double cxx[] = new double[order];
  double cy[]  = new double[order];
  double cyy[] = new double[order];
  double xmin = -3.0; // bounds for wind tunnel
  double xmax =  3.0;
  double ymin = -1.0;
  double ymax =  1.0;
  double xoff =  50.0;
  double yoff =  50.0;
  double scale = 100.0; // for plotting
  int nxg = 87; // 6.0 in X
  int nyg = 29; //  11 each side of 7 for cylinder, 2.0 in Y
  int ncyl = 37; // grid points in cylinder 2*3+2*5+3*7
  double hx = 6.0/86.0; // 87-1, x grid spacing
  double hy = 2.0/28.0; // 29-1, y grid spacing 
  double xg[] = new double[nxg]; // output grid
  double yg[] = new double[nyg];
  double u[][] = new double[nxg][nyg]; // velocity +X direction
  double v[][] = new double[nxg][nyg]; // velocity +Y direction
  double p[][] = new double[nxg][nyg]; // pressure
  int ixcyl[] = {40, 40, 40, 41, 41, 41, 41, 41, 42, 42, 42, 42, 42, 42, 42,
                 43, 43, 43, 43, 43, 43, 43, 44, 44, 44, 44, 44, 44, 44,
                 45, 45, 45, 45, 45, 46, 46, 46};
  int jxcyl[] = {13, 14, 15, 12, 13, 14, 15, 16, 11, 12, 13, 14, 15, 16, 17, 
                 11, 12, 13, 14, 15, 16, 17, 11, 12, 13, 14, 15, 16, 17,
                 12, 13, 14, 15, 16, 13, 14, 15}; 
  int nDOF = (nxg-2)*(nyg-2)-ncyl; // inside boundary, not in cylinder
             // number of degrees of freedom, u,v,p independent variable each
  int neqn = 3*nDOF; // number of equations 3*nDOF = number of variables
  int nlin = nDOF; //?? 3*order*order*nDOF; // maximum number of nonlinear terms
                            // uu for some uv, some vv
  int ntot= neqn+nlin; // max independent plus product variables neqn+nlin
  double A[][] = new double[neqn][ntot];
                    // A Y = X  system of non linear equations
  double Y[] = new double[neqn];
                    // RHS right hand side of equations
  double X[] = new double[ntot];
                    // initial guess and solution of non linear equations
  int var1[] = new int[ntot];    // variable number for non linear equations
  int var2[] = new int[ntot];    // variable number for non linear equations
  int var3[] = new int[ntot];    // variable number for non linear equations
  int iju[][] = new int[nxg][nyg];
  int ixpt[] = new int[nDOF]; // inverse of iju
  int jypt[] = new int[nDOF];

  int width = 700;  // window
  int height = 500; // upper part for title 
  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 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
  String fname;

// Non linear equations in A  nDOF u, v, p unknowns
// 0..nDOF-1            3nDOF-1   nlin       ntot = 3*nDOF+nlin
// |---u---|---v---|---p---|---uu, uv, vv---|
// |---var1[i]=i unknowns--|--var1, var2 ---|
//     var2, var3 = -1        variable index
// neqn rows ? may add continuity equations
// A*X=Y  Y is RHS
// function  incyl(i,j) = true for in cylinder, else false 
// function  bound(i,j) = 0 for not boundary, 1 for external, 2 for cylinder
// function  su(i,j) = subscript, index in A for u(x[i],y[j])  -1 for none
// function  sv(i,j) = subscript, index in A for v(x[i],y[j])  -1 for none
// function  sp(i,j) = subscript. index in A for p(x[i],y[j])  -1 for none
//

  cylinder_flow(String gname)
  {
    System.out.println("cylinder_flow.java running");
    fname = gname;
    System.out.println("using "+fname+", order="+order);
    build_vert();  

    setTitle(fname+" ");
    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());
    build_equations();

  } // end constructor cylinder_flow

  void build_vert()
  {
    double x, y;
    int ixnum = 0;

    // set up circle, 2D cylinder
    // xmin, xmax, ymin, ymax set manually for now
    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax);
    System.out.println("including boundary nxg="+nxg+
                       ", hx="+hx+", nyg="+nyg+", hy="+hy);

    // insert boundary and initial "guess"
    for(int i=0; i<nxg; i++) xg[i] = xmin+(double)i*hx; 
    System.out.println("xg[0]="+xg[0]+", xg[1]="+xg[1]+
                       ", xg["+(nxg-2)+"]="+xg[nxg-2]+
                       ", xg["+(nxg-1)+"]="+xg[nxg-1]);
    for(int j=0; j<nyg; j++) yg[j] = ymin+(double)j*hy; 
    System.out.println("yg[0]="+yg[0]+", yg[1]="+yg[1]+
                       ", yg["+(nyg-2)+"]="+yg[nyg-2]+
                       ", yg["+(nyg-1)+"]="+yg[nyg-1]);
    for(int i=0; i<nxg; i++)
    {
      for(int j=0; j<nyg; j++)
      {
	u[i][j] = u0;
        v[i][j] = v0;
        p[i][j] = p0;
        if(incyl(i,j))
	{
	  u[i][j] = 0.0;
	  v[i][j] = 0.0;
	  p[i][j] = 0.0;
	}
        if(bound(i,j)==0)
	{
	  iju[i][j] = ixnum; // for matrix A
          var1[ixnum] = ixnum; // for simeq_newton
          var1[ixnum+nDOF] = ixnum+nDOF;
          var1[ixnum+2*nDOF] = ixnum+2*nDOF;
          ixpt[ixnum] = i;
          jypt[ixnum] = j;
	  ixnum++;
	}
	else
	{
	  iju[i][j] = -1;
	}
        
      }
    }
    for(int i=0; i<ntot; i++)
    {
      var2[i] = -1;
      var3[i] = -1;
      if(i>neqn) var1[i] = -1;
    }
    System.out.println("check nDOF="+nDOF+"="+ixnum);
  } // end build_vert
    
  Boolean incyl(int i, int j) // return true if in airfoil cylinder
  {
    for(int k=0; k<ncyl; k++)
      if(i==ixcyl[k] && j==jxcyl[k]) return true;
    return false;
  } // end incyl

  int bound(int i, int j) // 0 for not boundary, 1 for external, 2 for airfoil
  {
    int jval;

    if(i==0 || i==nxg-1) return 1;
    if(j==0 || j==nyg-1) return 1;
    if(incyl(i,j)) return 2;
    return 0;
  } // end bound

  int su(int i, int j) // subscript, index in A for u(x[i],y[j])  -1 for none
  {
    return iju[i][j]; // built by  build_vert
  } // end su

  int sv(int i, int j) // subscript, index in A for v(x[i],y[j])  -1 for none
  {
    return iju[i][j]+nDOF;
  } // end sv

  int sp(int i, int j) // subscript, index in A for p(x[i],y[j])  -1 for none
  {
    return iju[i][j]+2*nDOF;
  } // end sp


  //  wind tunnel PDE, Navier Stokes static, incompressible, 2D
  //            area defined by xg[0]..xg[nxg-1]  yg[0]..yg[nyg-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 ?
  //            no z or t in this version
  //
  // 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 k=0;

    for(int i=1; i<nxg-1; i++)
    {
      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) cylinder airfoil edge
        x = xg[i];
        y = yg[j];
        new rderiv(1, order, k, hx, cx);
        new rderiv(2, order, k, hx, cxx);
        new rderiv(1, order, k, hy, cy);
        new rderiv(2, order, k, hy, cyy);

      }
    }

    // solve nonlinear equations, use results u, v, p
    new simeq_newton3(neqn, ntot-neqn, A, Y, var1, var2, var3, X,
                      6, 1.0, 1.0e-6, 0);
    //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 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,k;
    g.setColor(Color.black);
    g.drawString(fname, 20, 50);
    g.drawString("nxg="+nxg, 100, 50);
    g.drawString("nyg="+nyg, 200, 50);

    g.setColor(Color.red);
    for(int i=0; i<nxg; i++)
    {
      for(int j=0; j<nyg; j++)
      {
	k = bound(i,j);
        if(k==0) g.setColor(Color.yellow);
        if(k==1) g.setColor(Color.black);
        if(k==2) g.setColor(Color.blue);
 
        g.fillOval((int)(scale*(xg[i]-xmin)+xoff), 
	  	   height-(int)(scale*(yg[j]-ymin)+yoff),
                   3, 3);
      } // end j
    } // end i
  } // end paint

  public static void main(String args[])
  {
    String gname = "cylinder";
    // cylinder is just circle in two dimensions

    new cylinder_flow(gname);
    System.out.println("cylinder_flow.java ends");
  } // end main
} // end class cylinder_flow