/*
 * Couette flow - MAC(Marker and Cell) method 2D
 *
 * coded by Ikeuchi Mitsuru
 * ver 0.0.1   2006.02.25  created
 *
 */

import java.awt.*;
import java.awt.event.*;
import java.applet.*;

public class CouetteFlowMAC2D extends Applet
  implements ItemListener, ActionListener, AdjustmentListener, Runnable {

/* -------------------------------------- set global ------ */

  Choice ch_view;
  Button bt_reset, bt_startStop;
  Scrollbar sc_Reynolds,sc_vx0;
  Thread th = null;

  Dimension dim;
  Image imgOff;
  Graphics gOff;

  int sleepTime = 50;

  int resetSW = 0;
  int started = 1;
  int dispMode = 2;

  double t = 0.0;
  double dt = 0.01;
  double dx = 1.0/20.0;
  double dy = 1.0/20.0;
  double ReynoldsNumber = 100.0; 
  double alpha = 3.0;
  double vx0 = 1.0;
  double pressure0 = 1.0;
  double pError = 0.0;

  int NNx = 20;
  int NNy = 24;

  int cellKind[][] = new int[NNx][NNy];
  double vx[][] = new double[NNx][NNy];
  double vy[][] = new double[NNx][NNy];
  double nvx[][] = new double[NNx][NNy];
  double nvy[][] = new double[NNx][NNy];
  double source[][] = new double[NNx][NNy];
  double pressure[][] = new double[NNx][NNy];

/* ----------------------------- applet control ------ */

  public void init() {
    resize(630,360);
    setBackground(Color.white);

    dim = getSize();
    imgOff = createImage(dim.width,dim.height);
    gOff = imgOff.getGraphics();

    ch_view = new Choice();
    ch_view.addItem("flow");
    ch_view.addItem("pressure");
    ch_view.addItem("flow+P");
    ch_view.addItemListener(this);
    ch_view.select("flow+P");

    bt_reset= new Button("reset");
    bt_reset.addActionListener(this);

    bt_startStop= new Button("start/stop");
    bt_startStop.addActionListener(this); 
    
    sc_Reynolds = new Scrollbar(Scrollbar.HORIZONTAL,100,10,20,410);
    sc_Reynolds.addAdjustmentListener(this);

    sc_vx0 = new Scrollbar(Scrollbar.HORIZONTAL,100,10,0,110);
    sc_vx0.addAdjustmentListener(this);

    setLayout(new BorderLayout());
    Panel pnl = new Panel();
    pnl.setLayout(new GridLayout(1,6,5,0));
/*
    pnl.add(new Label("  "));
*/
    pnl.add(bt_reset);
    pnl.add(bt_startStop);
    pnl.add(sc_Reynolds);
    pnl.add(sc_vx0);
    pnl.add(new Label("  "));
    pnl.add(ch_view);
    add(pnl,"North");

    setInitialCondition();
  }

  public void start() {
    if (th == null) {
      th = new Thread(this);
      th.start();
    }
  }

  public void stop() {
    if (th != null) {
      th.stop();
      th = null;
    }
  }

  public void itemStateChanged(ItemEvent ev){
    if (ev.getSource() == ch_view) {
      dispMode = ch_view.getSelectedIndex();
    }
  }

  public void actionPerformed(ActionEvent ev){
    if(ev.getSource() == bt_reset){ 
      resetSW = 1;
    } else if (ev.getSource() == bt_startStop){
      if (started==0) {
        started = 1;
      } else {
        started = 0;
      }
    }
  }

  public void adjustmentValueChanged(AdjustmentEvent ev){
    if (ev.getSource() == sc_Reynolds) { 
      ReynoldsNumber = 1.0*(double)(sc_Reynolds.getValue());
    } else if(ev.getSource() == sc_vx0) { 
      vx0 = 0.01*(double)(sc_vx0.getValue());
	  setBoundary();
    }
  }

  public void run() {
    while (th != null) {
      try {
        timeEvolution();
        offPaint();
        repaint();
        Thread.sleep(sleepTime);
      } catch (InterruptedException e) { }
    }
  }

  public void update(Graphics g) {
    paint(g);
  }

  public void paint(Graphics g) {
    g.drawImage(imgOff,0,0,this);
  }

/* ----------------------------- offPaint -------------------- */

  void offPaint() {
    int gx,gy;
	
	gx = 100;
	gy = 80;

    gOff.setColor(Color.white);
    gOff.fillRect(0,0,dim.width,dim.height);

    if (dispMode==0) {
	  boundaryPlot(gx,gy);
      flowPlot(gx,gy, 200.0);

    } else if (dispMode==1) {
      pressurePlot(gx,gy, 50.0);
	  boundaryPlot(gx,gy);

    } else if (dispMode==2) {
      pressurePlot(gx,gy, 50.0);
	  boundaryPlot(gx,gy);
      flowPlot(gx,gy, 200.0);

    } else if (dispMode==3) {
      ;

    }

    gOff.setColor(Color.blue);
    gOff.drawString("t="+(int)(t*100.0+0.5)/100.0+" ",10,40);
    if (started==0) {
      gOff.drawString("stopped",630/6*1+10,40);
	} else {
      gOff.drawString("started",630/6*1+10,40);
	}
    gOff.drawString("Re="+(int)(ReynoldsNumber+0.5)+" ",630/6*2+10,40);
    gOff.drawString("vx0="+(int)(vx0*100.0+0.5)/100.0+" ",630/6*3+10,40);
	gOff.drawString("view",630/6*5+10,40);

    gOff.setColor(Color.black);
    gOff.drawString("Box="+(int)(NNx*dx*10.0+0.5)/10.0+"x"+(int)(NNy*dy*10.0+0.5)/10.0+" ",10,60);
/*
    gOff.setColor(Color.blue);
    gOff.drawString("dpMax="+(float)(pError)+" ",150,60);
*/
  }

/*------------------------- plot methods -------------------*/

  void flowPlot(int gx, int gy, double vmag) {
    int i,j,ix,iy,ix2,iy2;
    double sc,viewScale;

    sc = 10.0; viewScale=1.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        ix = (int)(sc*(i+0.5))+gx;
        iy = (int)(sc*(j+0.5))+gy;

        ix2 = ix+(int)(viewScale*vmag*vx[i][j]);
        iy2 = iy+(int)(viewScale*vmag*vy[i][j]);
	    if (vx[i][j]>0) {
	      gOff.setColor(Color.getHSBColor(0.70f,0.9f,0.8f));
	    } else {
          gOff.setColor(Color.getHSBColor(0.85f,0.9f,0.8f));
	    }
        gOff.drawLine(ix,iy, ix2,iy2);
      }
    }
  }

  void pressurePlot(int gx, int gy, double mag) {
    int i,j,ix,iy;
    double sc,viewScale,c,absc,p0;

    p0 = meanPressure();
    sc = 10.0; viewScale=1.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        ix = (int)(sc*i)+gx;
        iy = (int)(sc*j)+gy;
        c = mag*(pressure[i][j] - p0);
        if (c>0.999) c = 0.999;
        if (c<-0.999) c = -0.999;
        absc = c; if (absc<0) absc = -absc;
        gOff.setColor(Color.getHSBColor((float)(0.33-0.33*c),(float)(0.1+0.3*absc),0.8f));
        gOff.fillRect(ix,iy, 10,10);
      }
    }
  }

  void boundaryPlot(int gx, int gy) {
    int i,j,ix,iy;
    double sc;

    sc = 10.0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
	    if (cellKind[i][j]==1) {
          ix = (int)(sc*i)+gx;
          iy = (int)(sc*j)+gy;
	      gOff.setColor(Color.lightGray);
          gOff.fillRect(ix,iy, 10,10);
		}
      }
    }
  }

/*------------------------------------- statistics --------- */

 double meanPressure() {
    int i,j,n;
    double s;

    s = 0.0; n = 0;
    for (i=0; i<NNx; i++) {
      for (j=0; j<NNy; j++) {
        if (cellKind[i][j]==0) {
          s += pressure[i][j];
          n +=1;
        }
      }
    }
	return(s/n);
  }

/*-------------------------- set initial condition --------- */

  void setInitialCondition() {
    int i,j;
	
    t = 0.0;

    setBoundary();
    for (i=0; i<NNx; i++) {  
      for (j=0; j<NNy; j++) {
	    if (cellKind[i][j]==1) continue;
        pressure[i][j] = pressure0;
        vx[i][j] = 0.0;
        vy[i][j] = 0.0;
      }
    }
  }

  void setBoundary() {
    int i,j;

    for (i=0; i<NNx; i++) {  
      for (j=0; j<NNy; j++) {
	    if (j<2) {
	      cellKind[i][j] = 1;
          pressure[i][j] = pressure0;
          vx[i][j] = vx0;
          vy[i][j] = 0.0;
        } else if (j>=NNy-2) {
		  cellKind[i][j] = 1;
          pressure[i][j] = pressure0;
          vx[i][j] = 0.0;
          vy[i][j] = 0.0;
		} else {
          cellKind[i][j] = 0;
		}
      }
    }
  }

/* ----------------------------- timeEvolution -------------- */

  void timeEvolution() {

    if (resetSW==1) {
      setInitialCondition();
      resetSW = 0;
      started = 0;
    }

    if (started==1) {
      t = t + dt;
      vEvolution();
	  setNextPressure(30);
    }
  }


  void vEvolution() {
    int i,ip,ipp,im,imm,j,jp,jpp,jm,jmm;
    double absui,absvi, dpx,d2ux,d2uy,udux,vduy, dpy,d2vx,d2vy,udvx,vdvy;

    for (i=0; i<NNx; i++) {
	  ip = i+1; if (ip>=NNx) ip = 0;
	  ipp = ip+1; if (ipp>=NNx) ipp = 0;
	  im = i-1; if (im<0) im = NNx-1;
	  imm = im-1; if (imm<0) imm = NNx-1;	  
      for (j=0; j<NNy; j++) {
	    if (cellKind[i][j]==1) continue;
        jp = j+1; if (jp>=NNy) jp = 0;
        jpp = jp+1; if (jpp>=NNy) jpp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;
        jmm = jm-1; if (jmm<0) jmm = NNy-1;

        absui = vx[i][j]; if (absui<0.0) absui = -absui;
        absvi = vy[i][j]; if (absvi<0.0) absvi = -absvi;

        dpx = (pressure[ip][j]-pressure[im][j])/(2.0*dx);
        d2ux = (vx[ip][j]+vx[im][j]-2.0*vx[i][j])/(dx*dx);
        d2uy = (vx[i][jp]+vx[i][jm]-2.0*vx[i][j])/(dy*dy);

        udux = vx[i][j]*(-vx[ipp][j]+8.0*vx[ip][j]-8.0*vx[im][j]+vx[imm][j])/(12.0*dx) +
          alpha*absui*(vx[ipp][j]-4.0*vx[ip][j]+6*vx[i][j]-4.0*vx[im][j]+vx[imm][j])/(12.0*dx);
        vduy = vy[i][j]*(-vx[i][jpp]+8.0*vx[i][jp]-8.0*vx[i][jm]+vx[i][jmm])/(12.0*dy) +
          alpha*absvi*(vx[i][jpp]-4.0*vx[i][jp]+6*vx[i][j]-4.0*vx[i][jm]+vx[i][jmm])/(12.0*dy);
        nvx[i][j] = vx[i][j]+dt*(-udux-vduy-dpx+(d2ux+d2uy)/ReynoldsNumber);

		dpy = (pressure[i][jp]-pressure[i][jm])/(2.0*dy);
        d2vx = (vy[ip][j]+vy[im][j]-2.0*vy[i][j])/(dx*dx);
        d2vy = (vy[i][jp]+vy[i][jm]-2.0*vy[i][j])/(dy*dy);

        udvx = vx[i][j]*(-vy[ipp][j]+8.0*vy[ip][j]-8.0*vy[im][j]+vy[imm][j])/(12.0*dx) +
          alpha*absui*(vy[ipp][j]-4.0*vy[ip][j]+6*vy[i][j]-4.0*vy[im][j]+vy[imm][j])/(12.0*dx);
        vdvy = vy[i][j]*(-vy[i][jpp]+8.0*vy[i][jp]-8.0*vy[i][jm]+vy[i][jmm])/(12.0*dy) +
          alpha*absvi*(vy[i][jpp]-4.0*vy[i][jp]+6*vy[i][j]-4.0*vy[i][jm]+vy[i][jmm])/(12.0*dy);
        nvy[i][j] = vy[i][j]+dt*(-udvx-vdvy-dpy+(d2vx+d2vy)/ReynoldsNumber);
      }
    }

    for (i=0; i<NNx; i++) {  
      for (j=0; j<NNy; j++) {
	    if (cellKind[i][j]==1) continue;
        vx[i][j] = nvx[i][j];
        vy[i][j] = nvy[i][j];
      }
    }
  }

  void setNextPressure(int iter) {
    int ir;

    setSource();
	
    for (ir=0; ir<iter; ir++) {
      pError = poisson(1.5);
	}
  }

  void setSource() {
    int i,ip,im,j,jp,jm;
    double dux,dvx,duy,dvy;

    for (i=0; i<NNx; i++) {
	  ip = i+1; if (ip>=NNx) ip = 0;
	  im = i-1; if (im<0) im = NNx-1;	  
      for (j=0; j<NNy; j++) {
	    if (cellKind[i][j]==1) continue;
        jp = j+1; if (jp>=NNy) jp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;

        dux = (vx[ip][j]-vx[im][j])/(2.0*dx);
        dvx = (vy[ip][j]-vy[im][j])/(2.0*dx);
        duy = (vx[i][jp]-vx[i][jm])/(2.0*dy);
        dvy = (vy[i][jp]-vy[i][jm])/(2.0*dy);
		source[i][j] = -(dux+dvy)/dt +(dux*dux+2.0*dvx*duy+dvy*dvy);
      }
    }
  }

  double poisson(double omega) {
    int i,ip,im,j,jp,jm,ir,np;
    double dp,absdp,dpmax,p;

    setSource();

    dpmax = 0.0;
    for (i=0; i<NNx; i++) {
      ip = i+1; if (ip>=NNx) ip = 0;
      im = i-1; if (im<0) im = NNx-1;	  
      for (j=0; j<NNy; j++) {
        jp = j+1; if (jp>=NNy) jp = 0;
        jm = j-1; if (jm<0) jm = NNy-1;
        if (cellKind[i][j]==0) {
          dp = 0.25*(pressure[ip][j]+pressure[im][j]+pressure[i][jp]+pressure[i][jm]-4.0*pressure[i][j]+(dx*dx)*source[i][j]);
          pressure[i][j] += omega*dp;
          absdp = dp; if (absdp<0) absdp = -absdp;
          if (absdp>dpmax) dpmax = absdp;
        } else {
          p = 0.0; np = 0;
          if (cellKind[ip][j]==0) { p += pressure[ip][j]; np++; };
          if (cellKind[im][j]==0) { p += pressure[im][j]; np++; };
          if (cellKind[i][jp]==0) { p += pressure[i][jp]; np++; };
          if (cellKind[i][jm]==0) { p += pressure[i][jm]; np++; };
          if (np>0) {
            pressure[i][j] = p/np;
          } else {
            pressure[i][j] = pressure0;
          }
        }
      }
    }
    return ( dpmax );
  }


/* ----------------------------- end of applet -------------- */

}

/*
 * MAC(marker and cell) method 2D
 * 
 * imcompressible fluid
 * equation of continuity
 *   div(v)=0   (I)
 * conservation of momentum (Navier-Stokes)
 *   dv/dt + (v.nabra)v = -grad(p)+(1/Re)div(grad(v))   (II)
 *
 * finite difference
 * dv/dt --> {v(t+dt) - v(t)}/dt
 * dvx/dx --> {vx[i+1][j] - vx[i-1][j]}/(2dx),
 * dvx/dy --> {vx[i][j+1] - vx[i][j-1]}/(2dy),
 * ...
 * dvx2/dx2 --> {vx[i+1][j] + vx[i-1][j] - 2vx[i][j]}/(dx^2),
 *
 * Kawamura-Kuwahara scheme alpha = 3,
 *   u df/dx --> u[i][j]*dr + alpha*|u[i][j]|*ds,
 *     dr = {-f[i+2][j]+8f[i+1][j]-8f[i-1][j]+f[i-2][j]}/(12dx),
 *     ds = {f[i+2][j]-4f[i+1][j]+6f[i][j]-4f[i-1][j]+f[i-2][j]}/(12dx),
 *
 * pressure <-- (I),(II)
 * div(grad(p)) = dD/dt -{(dvx/dx)^2+(dvx/dy)(dvy/dx)+(dvy/dy)^2}
 *
 *
 * time evolution (vx,vy,p)
 *
 * (1) vx,vy,p -> next vx,vy
 *    v(t+dt) = v + dt(-v(v.nabra)v -grad(p)+(1/Re)div(grad(v)))
 *
 * (2) next vx,vy -> p(next)
 *    div(grad(p)) = -s,
 *      s = -dD/dt +{(dvx/dx)^2+(dvx/dy)(dvy/dx)+(dvy/dy)^2},
 *      D = div(v).
 *    solve Poissons equation 
 *
 * goto (1)
 *
 */