// cylinder_flow.c  cylinder in wind tunnel, uniform grid
//                  hard coded circle, two D view of cylinder
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "simeq_newton3.h"
#include "deriv.h"

static  char fname[] = "cylinder";
//static  const int order = 5; // 5th order approximation of derivatives
#define order 5
static  double cx[order][order];   // for rderiv
static  double cxx[order][order];  // first subscript is index
static  double cy[order][order];
static  double cyy[order][order];
static  double xmin = -2.0; // bounds for wind tunnel
static  double xmax =  2.0;
static  double ymin = -1.0;
static  double ymax =  1.0;
static  double xoff =  50.0;
static  double yoff =  50.0;
static  double scale = 100.0; // for plotting
// static  const int nxg = 61; // 4.0 in X
#define nxg 61
// static  const int nyg = 29; //  11 each side of 7 for cylinder, 2.0 in Y
#define nyg 29
// static  const int ncyl = 37; // grid points in cylinder 2*3+2*5+3*7
#define ncyl 37
static  double hx = 4.0/60.0; // 61-1, x grid spacing
static  double hy = 2.0/28.0; // 29-1, y grid spacing 
static  double xg[nxg]; // output grid
static  double yg[nyg];
static  double u[nxg][nyg]; // velocity +X direction
static  double v[nxg][nyg]; // velocity +Y direction
static  double p[nxg][nyg]; // pressure
static  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};
static  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}; 
// static  const int nDOF = (nxg-2)*(nyg-2)-ncyl; 
             // inside boundary, not in cylinder
             // number of degrees of freedom, u,v,p independent variable each
#define nDOF (nxg-2)*(nyg-2)-ncyl
// static  const int neqn = 3*nDOF; 
             // number of equations 3*nDOF = number of variables
#define neqn 3*nDOF
// static  const int nlin = 3*order*order*nDOF; guess ??? 
             // maximum number of nonlinear terms uu for some uv, some vv
#define nlin nDOF
// static  const int ntot = neqn+nlin; 
             // max independent plus product variables neqn+nlin
#define ntot neqn+nlin
static  double A[neqn][ntot];
                     // A Y = X  system of non linear equations
static  double Y[neqn];
                     // RHS right hand side of equations
static  double X[ntot];
                     // initial guess and solution of non linear equations
static  int var1[ntot];  // variable number for non linear equations
static  int var2[ntot];  // variable number for non linear equations
static  int var3[ntot];  // variable number for non linear equations
static  int iju[nxg][nyg];
static  int ixpt[nDOF]; // inverse of iju
static  int jypt[nDOF];

static  double u0 = 88.0;     // u velocity, uniform at input, output, top, bottom
static  double v0 = 0.0;      // v velocity is zero at these boundaries
static  double rho = 1.223;   // density
static  double mu = 182.0E-6; // viscosity
static  double p0 = 0.100;    // pressure at input and output, Pa about 14.5 psi

// 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
//

static void build_vert();
static int incyl(int i, int j); // return 1 if in airfoil cylinder
static int bound(int i, int j); // 0 for not boundary,
                                // 1 for external, 2 for airfoil
static int su(int i, int j);    // subscript, index in A for u(x[i],y[j])
                                //  -1 for none
static int sv(int i, int j);    // subscript, index in A for v(x[i],y[j])
                                //  -1 for none
static int sp(int i, int j);    // subscript, index in A for p(x[i],y[j])
                                // -1 for none
static void build_equations();

int main(int argc, char *argv[])
{
  int i;
  long long int alen;

  printf("cylinder_flow.c running \n");
  printf("%s  order=%d \n", fname, order);
  printf("including boundary nxg=%d, nyg=%d \n", nxg, nyg);
  printf("nDOF=%d, neqn=%d, nlin=%d, ntot=%d \n", nDOF, neqn, nlin, ntot);
  alen = (long long)8*(long long)neqn*(long long)ntot;
  printf("A = %lld bytes \n", alen);
  printf("xmin=%f, xmax=%f ymin=%f, ymax=%f \n",
	 xmin, xmax, ymin, ymax);
  printf("hx=%f, hy=%f \n", hx, hy);
  build_vert();  
  build_equations();
  printf("starting to solve non-linear system of equations %d by %d \n",
         neqn, ntot-neqn);
  fflush(stdout);
  // solve nonlinear equations, use results u, v, p
  //?? simeq_newton3(neqn, ntot-neqn, A, Y, var1, var2, var3, X,
  //??                    6, 1.0e-6, 0);
  // print results from X
  printf(" ix  jy    u            v            p \n");
  for(i=0; i<nDOF; i++)
  {
    printf("%3d %3d %e %e %e \n", ixpt[i], jypt[i],
           X[i], X[i+nDOF], X[i+2*nDOF]);
  }

  printf("cylinder_flow.c finished \n");
  return 0;
} // end main


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

  // set up circle, 2D cylinder
  // xmin, xmax, ymin, ymax set manually for now
  // insert boundary and initial "guess"
  for(i=0; i<nxg; i++) xg[i] = xmin+(double)i*hx; 
  printf("xg[0]=%f, xg[1]=%f, xg[%2d]=%f, xg[%2d]=%f \n",
         xg[0], xg[1], nxg-2, xg[nxg-2], nxg-1, xg[nxg-1]);
  for(j=0; j<nyg; j++) yg[j] = ymin+(double)j*hy; 
  printf("yg[0]=%f, yg[1]=%f, yg[%2d]=%f, yg[%2d]=%f \n",
         yg[0], yg[1], nyg-2, yg[nyg-2], nyg-1, yg[nyg-1]);
  for(i=0; i<nxg; i++)
  {
    for(j=0; j<nyg; j++)
    {
      u[i][j] = u0;
      v[i][j] = v0;
      p[i][j] = p0;
      if(incyl(i,j)!=0)
      {
        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(i=0; i<ntot; i++)
  {
    var2[i] = -1;
    var3[i] = -1;
    if(i>neqn) var1[i] = -1;
  }
  printf("check nDOF=%d = %d ixnum \n", nDOF, ixnum);
} // end build_vert
    
static int incyl(int i, int j) // return 1 if in airfoil cylinder
{
  int k;

  for(k=0; k<ncyl; k++)
    if(i==ixcyl[k] && j==jxcyl[k]) return 1;
  return 0;
} // end incyl

static int bound(int i, int j) // 0 for not boundary,
                               // 1 for external, 2 for airfoil
{
  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

static 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

static 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

static 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 

static void build_equations()
{
  double x, y;
  int k=0;
  int i, j, ix, jy, iua, iva, ipa;
  int urow, vrow, prow;

  printf("derivatives \n");
  for(k=0; k<order; k++)
  {
    rderiv(1, order, k, hx, cx[k]);
    rderiv(2, order, k, hx, cxx[k]);
    rderiv(1, order, k, hy, cy[k]);
    rderiv(2, order, k, hy, cyy[k]);
    for(i=0; i<order; i++)
      printf("k=%d, i=%d, cx=%f, cxx=%f, cy=%f, cyy=%f \n",
             k, i, cx[k][i],  cxx[k][i],  cy[k][i],  cyy[k][i]); 
    printf(" \n");
  }
  printf(" \n");

  for(i=1; i<nxg-1; i++)
  {
    for(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];
      iua = su(i,j);
      iva = sv(i,j);
      ipa = sp(i,j);
      ix = ixpt[iua];
      jy = jypt[iua];
      printf("i=%d, j=%d, x=%f, y=%f, ix=%d, jy=%d, iua, iva, ipa \n",
             i, j, x, y, ix, jy, iua, iva, ipa);

      // build A and Y=RHS, check for near boundary
      urow = su(i,j); // also column for u
      vrow = sv(i,j); // also column for v
      prow = sp(i,j); // also column for p
      // special cases near boundary X direction
      if(i==1) // cx[0][]  cxx[0][]
      {
      }
      else if(i==nxg-2) // cx[order-1][]
      {
      }
      else if(bound(i-1,j)!=0) // cx[0][]
      {
      }
      else if(bound(i+1,j)!=0) // cx[order-1][]
      {
      }
      else // cx[2] [(order+1)/2]
      {
        for(k=0; k<order; k++)
        {
          A[urow][urow+k] = cx[2][k]; // ?? need actual equation, cxx
          A[vrow][vrow+k] = cy[2][k]; // ?? need actual equation
          A[prow][prow+k] = cx[2][k]; // ?? need actual equation
        }
      }
      // special cases near boundary Y direction
      if(j==1) // cy[0][]  cy[00][]
      {
      }
      else if(j==nyg-2) // cy[order-1][]
      {
      }
      else if(bound(i,j-1)!=0) // cy[0][]
      {
      }
      else if(bound(i,j+1)!=0) // cy[order-1][]
      {
      }
      else // cy[2] cy[(order+1)/2]
      {
        for(k=0; k<order; k++)
        {
          A[vrow][vrow+k] = cy[2][k]; // ?? need actual equation
          A[prow][prow+k] = cx[2][k]; // ?? need actual equation
        }
      }
    }
  }
}  // end build_equations

// end cylinder_flow.c