//Natural Convection 
//Android Version of the Lattice Boltzmann Method
//Copyright @2014- Takeshi Seta All Rights Reserved.
package jp.seta.lbm_natural_convection;

import android.os.Bundle;
import android.app.Activity;
import android.view.View;
import android.graphics.*;
import android.content.Context;
import android.os.Handler;
import android.os.Message;
import android.view.MotionEvent;

public class MainActivity extends Activity {
   LBMDrawView view;
 
    @Override
	public void onCreate(Bundle savedInstanceState) {
        super.onCreate(savedInstanceState);
        LBMDrawView v = new LBMDrawView(getApplication());
        setContentView(v);
	}
}

//Animation
class RedrawHandler extends Handler {
    private View view;
    private int delayTime;
    private int frameRate;
    public RedrawHandler(View view, int frameRate) {
        this.view = view;
        this.frameRate = frameRate;
    }
    public void start() {
        this.delayTime = 1000 / frameRate;
        this.sendMessageDelayed(obtainMessage(0), delayTime);
    }
    public void stop() {
        delayTime = 0;
    }
    @Override
    public void handleMessage(Message msg) {
        view.invalidate();
        if (delayTime == 0) return; // stop
        sendMessageDelayed(obtainMessage(0), delayTime);
    }
}

class LBMDrawView extends View {
    public int nx =  51, ny =  51, time = 0, in, jn, tx = 280, ty = 100;
    public double  wid = 8d, vc = 0.2d, beta, pr, ra, gx, gy, kap, nu, umax, vmax;
    public double tauf = 0.8d, taug = 0.8d, el = 1d, er = 0d, et, eb, tmp, u2, norm;
    double[][] u = new double[nx][ny], v = new double[nx][ny], rho = new double[nx][ny];
    double[][] un= new double[nx][ny], vn= new double[nx][ny];
    double[][] x = new double[nx][ny], y = new double[nx][ny];
    double[][][] f = new double[9][nx][ny], f0 = new double[9][nx][ny], ftmp = new double[9][nx][ny];
    double[] cx = new double[9], cy = new double[9], nu2 = new double[nx];
    double[][] e = new double[nx][ny], en= new double[nx][ny], q= new double[nx][ny];
    double[][][] g = new double[9][nx][ny], g0 = new double[9][nx][ny], gtmp = new double[9][nx][ny], fi = new double[9][nx][ny];

    Paint p = new Paint();
    Paint t = new Paint();
        
    public boolean onTouchEvent(MotionEvent event) {
        tx = (int)event.getX();
        ty = (int)event.getY();
        if(tx >= 430){ tx = 430;}
        if(tx <= 130){ tx = 130;}
        return true;
    }    
    
    public LBMDrawView(Context c) {
        super(c);
        // initial condition
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
    	    u[i][j] = 0.0d; v[i][j] = 0.0d; rho[i][j] = 1.0d; 
    	    x[i][j] = i*wid + 50; y[i][j] = j*wid + 170;
            e[i][j] = (double)(nx - 1 - i)/(double)(nx - 1);
        } }
        
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            un[i][j] = u[i][j]; vn[i][j] = v[i][j]; en[i][j] = e[i][j];
        } }
        
        gx = 0d; gy = 0d;
        ra = 10000d; pr = 0.71d;

        beta = vc*vc/(double)(nx-1);
        tauf = 3d   *vc*(double)(nx-1)*Math.sqrt(pr/ra) + 0.5d;
        taug = 3d/2d*vc*(double)(nx-1)/Math.sqrt(pr*ra) + 0.5d;
        kap = 2d/3d*(taug - 0.5d);

        cx[0] = 0d; cy[0] = 0d;
        cx[1] = 1d; cy[1] = 0d; cx[2] = 0d; cy[2] = 1d;
        cx[3] =-1d; cy[3] = 0d; cx[4] = 0d; cy[4] =-1d;
        cx[5] = 1d; cy[5] = 1d; cx[6] =-1d; cy[6] = 1d;
        cx[7] =-1d; cy[7] =-1d; cx[8] = 1d; cy[8] =-1d;

        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            u2 = u[i][j]*u[i][j] + v[i][j]*v[i][j];      
            f[0][i][j] = rho[i][j]*(1d - 3d/2d*u2)*4d/9d;
            g[0][i][j] =-rho[i][j]*e[i][j]*u2*2d/3d;
            for (int k = 1; k < 5; k++) {
                tmp = cx[k]*u[i][j] + cy[k]*v[i][j];      
                f[k][i][j] = rho[i][j]*(1d + 3d*tmp + 9d/2d*tmp*tmp - 3d/2d*u2)/9d;
                g[k][i][j] = rho[i][j]*e[i][j]*(3d/2d + 3d/2d*tmp + 9d/2d*tmp*tmp - 3d/2d*u2)/9d;
            }
            for (int k = 5; k < 9; k++) {
                tmp = cx[k]*u[i][j] + cy[k]*v[i][j];
                f[k][i][j] = rho[i][j]*(1d + 3d*tmp + 9d/2d*tmp*tmp - 3d/2d*u2)/36d;
                g[k][i][j] = rho[i][j]*e[i][j]*(3d + 6d*tmp + 9d/2d*tmp*tmp - 3d/2d*u2)/36d;
            }
        } }
        setFocusable(true);
        RedrawHandler handler = new RedrawHandler(this, 8);
        handler.start();        
    }
    
    public void onDraw(Canvas c) { 
        ra = Math.pow(10d,((double)tx - 130d)/(430d - 130d)*2d + 3d);
        tauf = 3d   *vc*(double)(nx-1)*Math.sqrt(pr/ra) + 0.5d;
        taug = 3d/2d*vc*(double)(nx-1)/Math.sqrt(pr*ra) + 0.5d;
        kap = 2d/3d*(taug - 0.5d);

        norm = 0.0d;
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            tmp = Math.abs(e[i][j] - en[i][j]);
            if(tmp > norm) norm = tmp;
            tmp = Math.abs(Math.sqrt( u[i][j]* u[i][j] +  v[i][j]* v[i][j])
                          -Math.sqrt(un[i][j]*un[i][j] + vn[i][j]*vn[i][j]));
            if(tmp > norm) norm = tmp;
        } }
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            un[i][j] = u[i][j]; vn[i][j] = v[i][j]; en[i][j] = e[i][j];
        } }

    	//forcing
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            gx = 0d; gy = - beta*(e[i][j] - 0.5d);
            fi[0][i][j] = 0d;
            for (int k = 1; k < 5; k++) {
              fi[k][i][j] = (gx*cx[k] + gy*cy[k])/ 3d*rho[i][j];
            }
            for (int k = 5; k < 9; k++) {
              fi[k][i][j] = (gx*cx[k] + gy*cy[k])/12d*rho[i][j];
            }
        } }

    	//collision
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            u2 = u[i][j]*u[i][j] + v[i][j]*v[i][j];      
            f0[0][i][j] = rho[i][j]*(1d - 3d/2d*u2)*4d/9d;
            g0[0][i][j] =-rho[i][j]*e[i][j]*u2*2d/3d;
            for (int k = 1; k < 5; k++) {
                tmp = cx[k]*u[i][j] + cy[k]*v[i][j];      
                f0[k][i][j] = rho[i][j]*(1d + 3d*tmp + 9d/2d*tmp*tmp - 3d/2d*u2)/9d;
                g0[k][i][j] = rho[i][j]*e[i][j]*(3d/2d + 3d/2d*tmp + 9d/2d*tmp*tmp - 3d/2d*u2)/9d;
            }
            for (int k = 5; k < 9; k++) {
                tmp = cx[k]*u[i][j] + cy[k]*v[i][j];
                f0[k][i][j] = rho[i][j]*(1d + 3d*tmp + 9d/2d*tmp*tmp - 3d/2d*u2)/36d;
                g0[k][i][j] = rho[i][j]*e[i][j]*(3d + 6d*tmp + 9d/2d*tmp*tmp - 3d/2d*u2)/36d;
            }
        } }
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) { for(int k = 0; k < 9; k++){
            f[k][i][j] = f[k][i][j] - (f[k][i][j] - f0[k][i][j])/tauf + fi[k][i][j];
            g[k][i][j] = g[k][i][j] - (g[k][i][j] - g0[k][i][j])/taug;
        } } }
            
      	// streaming
    	for (int i = 0; i < nx; i++){ for (int j = 0; j < ny; j++){ for (int k = 0; k < 9; k++){
    	    ftmp[k][i][j] = f[k][i][j];
    	    gtmp[k][i][j] = g[k][i][j];
        } } }
        for (int i = 0; i < nx-1; i++){ for (int j = 0; j < ny; j++){
            in = i + 1;
            f[1][in][j] = ftmp[1][i][j];
            g[1][in][j] = gtmp[1][i][j];
        } }
        for (int i = 0; i < nx; i++){ for (int j = 0; j < ny-1; j++){
            jn = j + 1;
            f[2][i][jn] = ftmp[2][i][j];
            g[2][i][jn] = gtmp[2][i][j];
        } }
        for (int i = 1; i < nx; i++){ for (int j = 0; j < ny; j++){
            in = i - 1;
            f[3][in][j] = ftmp[3][i][j];
            g[3][in][j] = gtmp[3][i][j];
        } }
        for (int i = 0; i < nx; i++){ for (int j = 1; j < ny; j++){
            jn = j - 1;
            f[4][i][jn] = ftmp[4][i][j];
            g[4][i][jn] = gtmp[4][i][j];
        } }
        for (int i = 0; i < nx-1; i++){ for (int j = 0; j < ny-1; j++){
            in = i + 1;  jn = j + 1;
            f[5][in][jn] = ftmp[5][i][j];
            g[5][in][jn] = gtmp[5][i][j];
        } }
        for (int i = 1; i < nx; i++){ for (int j = 0; j < ny-1; j++){
            in = i - 1; jn = j + 1;
            f[6][in][jn] = ftmp[6][i][j];
            g[6][in][jn] = gtmp[6][i][j];
        } }
        for (int i = 1; i < nx; i++){ for (int j = 1; j < ny; j++){
            in = i - 1; jn = j - 1;
            f[7][in][jn] = ftmp[7][i][j];
            g[7][in][jn] = gtmp[7][i][j];
        } }
        for (int i = 0; i < nx-1; i++){ for (int j = 1; j < ny; j++){
            in = i + 1;  jn = j - 1;
            f[8][in][jn] = ftmp[8][i][j];
            g[8][in][jn] = gtmp[8][i][j];
        } }

        //boundary condition
        // Q. Zou and X. He, Phys. Fluids 9, 1591 (1997).
        for (int j = 0; j < ny; j++){
      	    rho[0][j] =f[0][0][j] + f[2][0][j] + f[4][0][j]
                + 2d*(f[6][0][j] + f[3][0][j] + f[7][0][j]);
            f[1][0][j] = f[3][0][j];
            f[5][0][j] = f[7][0][j] + 0.5d*(f[4][0][j] - f[2][0][j]);
            f[8][0][j] = f[6][0][j] - 0.5d*(f[4][0][j] - f[2][0][j]);
            g[1][0][j] =-g[3][0][j] +      (f[1][0][j] + f[3][0][j]) + rho[0][j]*(el/3d - 2d/9d);
            g[5][0][j] =-g[7][0][j] + 2.0d*(f[5][0][j] + f[7][0][j]) + rho[0][j]*(el/6d - 1d/9d);
            g[8][0][j] =-g[6][0][j] + 2.0d*(f[8][0][j] + f[6][0][j]) + rho[0][j]*(el/6d - 1d/9d);
            
            rho[nx-1][j] =f[0][nx-1][j] + f[2][nx-1][j] + f[4][nx-1][j]
                    + 2d*(f[5][nx-1][j] + f[1][nx-1][j] + f[8][nx-1][j]);
            f[3][nx-1][j] = f[1][nx-1][j];
            f[7][nx-1][j] = f[5][nx-1][j] + 0.5d*(f[2][nx-1][j] - f[4][nx-1][j]);
            f[6][nx-1][j] = f[8][nx-1][j] - 0.5d*(f[2][nx-1][j] - f[4][nx-1][j]);
            g[3][nx-1][j] =-g[1][nx-1][j] +      (f[3][nx-1][j] + f[1][nx-1][j]) + rho[nx-1][j]*(er/3d - 2d/9d);
            g[7][nx-1][j] =-g[5][nx-1][j] + 2.0d*(f[7][nx-1][j] + f[5][nx-1][j]) + rho[nx-1][j]*(er/6d - 1d/9d);
            g[6][nx-1][j] =-g[8][nx-1][j] + 2.0d*(f[6][nx-1][j] + f[8][nx-1][j]) + rho[nx-1][j]*(er/6d - 1d/9d);
        }

        for (int i = 0; i < nx; i++){
            rho[i][ny-1] =f[0][i][ny-1] + f[1][i][ny-1] + f[3][i][ny-1]
                    + 2d*(f[2][i][ny-1] + f[5][i][ny-1] + f[6][i][ny-1]);
            f[4][i][ny-1] = f[2][i][ny-1];
            f[7][i][ny-1] = f[5][i][ny-1] + 0.5d*(f[1][i][ny-1] - f[3][i][ny-1]);
            f[8][i][ny-1] = f[6][i][ny-1] - 0.5d*(f[1][i][ny-1] - f[3][i][ny-1]);

            et = 4d*e[i][ny-2]/3d - e[i][ny-3]/3d; 
            g[4][i][ny-1] =-g[2][i][ny-1] +      (f[4][i][ny-1] + f[2][i][ny-1]) + rho[i][ny-1]*(et/3d - 2d/9d);
            g[7][i][ny-1] =-g[5][i][ny-1] + 2.0d*(f[7][i][ny-1] + f[5][i][ny-1]) + rho[i][ny-1]*(et/6d - 1d/9d);
            g[8][i][ny-1] =-g[6][i][ny-1] + 2.0d*(f[8][i][ny-1] + f[6][i][ny-1]) + rho[i][ny-1]*(et/6d - 1d/9d);
            
            rho[i][0] = f[0][i][0] + f[1][i][0] + f[3][i][0]
                  + 2d*(f[4][i][0] + f[7][i][0] + f[8][i][0]);
            f[2][i][0] = f[4][i][0];
            f[5][i][0] = f[7][i][0] - 0.5d*(f[1][i][0] - f[3][i][0]);
            f[6][i][0] = f[8][i][0] + 0.5d*(f[1][i][0] - f[3][i][0]);

            eb = 4d*e[i][1]/3d - e[i][2]/3d; 
	        g[2][i][0] =-g[4][i][0] +      (f[2][i][0] + f[4][i][0]) + rho[i][0]*(eb/3d - 2d/9d);
            g[5][i][0] =-g[7][i][0] + 2.0d*(f[5][i][0] + f[7][i][0]) + rho[i][0]*(eb/6d - 1d/9d);
            g[6][i][0] =-g[8][i][0] + 2.0d*(f[7][i][0] + f[8][i][0]) + rho[i][0]*(eb/6d - 1d/9d);
        }

        //physics
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            rho[i][j] = f[0][i][j]; u[i][j] = 0.0d; v[i][j] = 0.0d; e[i][j] = g[0][i][j];
            for (int k = 1; k < 9; k++){
                rho[i][j] = rho[i][j] + f[k][i][j];
                  u[i][j] =   u[i][j] + f[k][i][j] * cx[k];
                  v[i][j] =   v[i][j] + f[k][i][j] * cy[k];
                  e[i][j] =   e[i][j] + g[k][i][j];
            }
            if(rho[i][j] != 0){
                u[i][j] = u[i][j]/rho[i][j];
                v[i][j] = v[i][j]/rho[i][j];
                e[i][j] = e[i][j]/rho[i][j];
            }else{
                u[i][j] = 0d; v[i][j] = 0d; e[i][j] = 0d;
            }
        } }
        
        //maximum velocity
      	umax = 0d; vmax = 0d;
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            if(Math.abs(u[i][j]) > umax){ umax = Math.abs(u[i][j]); }
            if(Math.abs(v[i][j]) > vmax){ vmax = Math.abs(v[i][j]); }
        } }
      	umax = umax*(double)(nx-1)/kap;
      	vmax = vmax*(double)(nx-1)/kap;
        
        //Nusselt number
        for (int j = 0; j < ny; j++) {
 	      q[   0][j] = (-e[   2][j] + 4d*e[   1][j] - 3d*e[   0][j])/2d;
 	      q[nx-1][j] = ( e[nx-3][j] - 4d*e[nx-2][j] + 3d*e[nx-1][j])/2d;
          for (int i = 1; i < nx-1; i++) {
            q[i][j] = (e[i+1][j] - e[i-1][j])/2d;
          }
        }
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
 	        q[i][j] = u[i][j]*e[i][j] - kap*q[i][j];
        } }

        for (int i = 0; i < nx; i++) {
          nu2[i] = 0.0d;
        }
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny-2; j = j+2) {
          nu2[i] = nu2[i] + (q[i][j] + 4d*q[i][j+1] + q[i][j+2])/3d;
        } }

        nu = 0.0d;
        for (int i = 0; i < nx-2; i = i+2) {
          nu = nu + (nu2[i] + 4d*nu2[i+1] + nu2[i+2])/3d;
        }
        nu = nu/kap/(double)(nx - 1);

        //graphics
      	c.drawColor(Color.WHITE);
      	tmp = 0d;
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            u2 = (float)Math.sqrt(u[i][j]*u[i][j] + v[i][j]*v[i][j]);      
            if(u2 > tmp){ tmp = u2; }
        } }

      	p.setStrokeWidth(3);
        p.setStyle(Paint.Style.STROKE);
        p.setColor(Color.GRAY);
        c.drawRect((float)(x[0][0]), (float)(y[0][0]), (float)(x[nx-1][ny-1]+ wid), (float)(y[nx-1][ny-1]+ wid), p);            
   	    
 	    
        p.setStrokeWidth(2);
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
          // temperature
          in = (int)(e[i][j]*255); jn =  (int)((1-e[i][j])*255);
      	  p.setColor(Color.argb(127,in,(int)(e[i][j]*127),jn));
          p.setStyle(Paint.Style.FILL);
          c.drawRect((float)(x[i][j]), (float)(y[i][j]), (float)(x[i][j] + wid), (float)(y[i][j] + wid), p);            

       	  // velocity
          p.setColor(Color.argb(255,in,(int)(e[i][j]*127),jn));
          c.drawLine((float)(x[i][j] + wid/2)                     , (float)(y[i][j] + wid/2),
                     (float)(x[i][j] + wid/2 + u[i][j]/tmp*wid*10), (float)(y[i][j] + wid/2 + v[i][j]/tmp*wid*10), p); 
        } }

        p.setColor(Color.argb(255,0, 255, 255));
      	p.setStrokeWidth(3);
        c.drawLine(130f, 120f, 430f, 120f, p); 
        c.drawLine(130f, 115f, 130f, 130f, p); 
        c.drawLine(280f, 115f, 280f, 130f, p); 
        c.drawLine(430f, 115f, 430f, 130f, p); 

      	p.setStrokeWidth(1);
        p.setColor(Color.RED);
        c.drawLine((float)(tx)-4, 140f, (float)(tx)-4, 124f,p);
        c.drawLine((float)(tx)-3, 140f, (float)(tx)-3, 123f,p);
        c.drawLine((float)(tx)-2, 140f, (float)(tx)-2, 122f,p);
        c.drawLine((float)(tx)-1, 140f, (float)(tx)-1, 121f,p);
        c.drawLine((float)(tx)  , 140f, (float)(tx)  , 120f,p);
        c.drawLine((float)(tx)+1, 140f, (float)(tx)+1, 121f,p);
        c.drawLine((float)(tx)+2, 140f, (float)(tx)+2, 122f,p);
        c.drawLine((float)(tx)+3, 140f, (float)(tx)+3, 123f,p);
        c.drawLine((float)(tx)+4, 140f, (float)(tx)+4, 124f,p);
        
        // Text
        t.setColor(Color.BLACK);
        t.setTextSize (16);
        String msg1 = "Natural"+" convection"+", 51x51" +" grid" + ", Pr = 0.71";
        c.drawText(msg1, 10, 30, t);
        String msg2 = "Time: " +time++ +", Norm: " +norm;
        c.drawText(msg2, 10, 50, t);

        t.setColor(Color.BLACK);
        String msg3 = "Nusselt number: "+nu;
        c.drawText(msg3, 10, 70, t);

        t.setColor(Color.BLACK);
        String msg4 = "Maximum Velocity: ("+(float)umax+", "+(float)vmax+")";
        c.drawText(msg4, 10, 90, t);
        
        t.setColor(Color.BLUE);
        String msg5 =   "1000"; c.drawText(msg5, 110, 110, t);
        String msg6 =  "10000"; c.drawText(msg6, 260, 110, t);
        String msg7 = "100000"; c.drawText(msg7, 400, 110, t);
        
        t.setColor(Color.BLACK);
        String msg8 = "Ra:";
        c.drawText(msg8, 10, 125, t);

        t.setColor(Color.RED);
        String msg9 = ""+(float)ra;
        c.drawText(msg9, 40, 125, t);

        t.setColor(Color.BLUE);
        String msg10 = "Change the Rayleigh number.";
        c.drawText(msg10, 170, 157, t);

        t.setColor(Color.BLACK);
        String msg11 = "Y."+" Peng,"+" C."+" Shu,"+" and"+" Y."+" T."+" Chew,"+" Phys."+" Rev."+" E,"+" 68,"+" 026701"+" (2003).";
        c.drawText(msg11, 10, 600, t);
        String msg12 = "Android Version of the Lattice Boltzmann Method";
        c.drawText(msg12, 10, 620, t);
        String msg13 = "Copyright @2014- Takeshi Seta All Rights Reserved.";
        c.drawText(msg13, 10, 640, t);
    }
}