// Particle motion induced by a circular ripple 
// Android Version of the Lattice Boltzmann Method combined with smoothed-profile method
// Copyright @2014- Takeshi Seta, Kento Tanaka, Masahiro Matsuzuru, Takayuki Yamagishi All Rights Reserved.
       
package jp.seta.LBM_Particle;

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 = 10 / 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 = 32, ny = 32, in, jn, tx, ty, tflag = 0, np = 5;
    public double tau = 1.7d, mu, tmp, u2, wid = 14.0d, u0 = 0.1d, dt, time = 0.0d;
    double[][] u = new double[nx][ny], v = new double[nx][ny], rho = new double[nx][ny];
    double[][] x = new double[nx][ny], y = new double[nx][ny];
    double[][] fx = new double[nx][ny], fy = new double[nx][ny];
    double[][] up = new double[nx][ny], vp = new double[nx][ny], phi = new double[nx][ny];
    double[][][] f = new double[9][nx][ny], f0 = new double[9][nx][ny];
    double[][][]  ftmp = new double[9][nx][ny], fi = new double[9][nx][ny];
    double[] cx = new double[9], cy = new double[9];
    
    public double tmp1, tmp2, rhop, guzai, mp, iip, rp, zetaw, zetap, gxp, gyp;
    double[] theta = new double[np], fxp = new double[np], fyp = new double[np];
    double[] xp = new double[np], yp = new double[np], fxc = new double[np], fyc = new double[np];
    double[] omega = new double[np],  omega1 = new double[np], omega2  = new double[np], torqp = new double[np];
    double[] upt = new double[np], up1 = new double[np], up2 = new double[np];
    double[] vpt = new double[np], vp1 = new double[np], vp2 = new double[np];
    		
    Paint p = new Paint();
    Paint t = new Paint();

    public boolean onTouchEvent(MotionEvent event) {
        tx = (int)event.getX();
        ty = (int)event.getY();
        tx = (int)(tx/wid);
        ty = (int)(ty/wid);
        tflag = 1;
        return true;
    }

    public LBMDrawView(Context c) {
        super(c);
        // initial condition

        mu = (tau - 0.5d)/3.0d;
        dt = mu/((double)(nx)/0.5d)/((double)(nx)/0.5d)/0.1d;
        
        rp  = (double)nx*0.030625d/0.5d;
        rhop = 0.7d;
        mp = Math.PI*rp*rp*rhop;
        iip = Math.PI*rp*rp*rp*rp*0.5d*rhop;
        guzai = 2.0d;
        zetaw = 1.0d;
        zetap = 1.5d;
        gxp = 0.0d;
        gyp = 981d/0.5d*(double)(nx)*dt*dt;
        gyp = 0.0d;

        xp[0] = (double)nx/2.0d  ; yp[0] = (double)ny/2.0d;
        xp[1] = xp[0] + 2*rp + 4d; yp[1] = yp[0];
        xp[2] = xp[0] - 2*rp - 4d; yp[2] = yp[0];
        xp[3] = xp[0]            ; yp[3] = yp[0] - 2*rp - 4d;
        xp[4] = xp[0]            ; yp[4] = yp[0] + 2*rp + 4d;
        
        for (int m = 0; m < np; m++) {
            upt[m] = 0.0d;    up1[m] = 0.0d;    up2[m] = 0.0d;
            vpt[m] = 0.0d;    vp1[m] = 0.0d;    vp2[m] = 0.0d;
          omega[m] = 0.0d; omega1[m] = 0.0d; omega2[m] = 0.0d;
          theta[m] = 0.0d;  torqp[m] = 0.0d;
          fxp[m] = 0.0d;      fyp[m] = 0.0d;
          fxc[m] = 0.0d;      fyc[m] = 0.0d;
        }

        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] = (double)(i)*wid; y[i][j] = (double)(j)*wid;
        } }

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

        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]*(1.0d - 3.0d/2.0d*u2)*4.0d/9.0d;
            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]*(1.0d + 3.0d*tmp + 9.0d/2.0d*tmp*tmp - 3.0d/2.0d*u2)/9.0d;
            }
            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]*(1.0d + 3.0d*tmp + 9.0d/2.0d*tmp*tmp - 3.0d/2.0d*u2)/36.0d;
            }
         } }

        //S-function
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++){
            up[i][j] = 0.0d; vp[i][j] = 0.0d; phi[i][j] = 0.0d;
            for (int m = 0; m < np; m++) {
              tmp1 = rp - Math.sqrt((xp[m] - (double)(i))*(xp[m] - (double)(i)) + (yp[m] - (double)(j))*(yp[m] - (double)(j)));
              if(tmp1 > guzai/2.0d){
                tmp2 = 1.0d;
              }else if(tmp1 < - guzai/2.0d){
                tmp2 = 0.0d;
              }else{
                tmp2 = (Math.sin(Math.PI*tmp1/guzai) + 1.0d)/2.0d;
              }
              up[i][j] =  up[i][j] + tmp2*(upt[m] - omega[m]*((double)(j) - yp[m]));
              vp[i][j] =  vp[i][j] + tmp2*(vpt[m] + omega[m]*((double)(i) - xp[m]));
   	         phi[i][j] = phi[i][j] + tmp2;
            }
        } }

        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++){
            if(phi[i][j] != 0){
                up[i][j] = up[i][j]/phi[i][j];
                vp[i][j] = vp[i][j]/phi[i][j];
   	     }else{
                up[i][j] = 0.0d;
                vp[i][j] = 0.0d;
             }
        } }

        setFocusable(true);
        RedrawHandler handler = new RedrawHandler(this, 8);
        handler.start();
    }

    public void onDraw(Canvas c) {
	//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)*4f/9d;
            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 + 3f*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;
            }
        } }

        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++)  { for (int k = 0; k < 9; k++) {
   	            fi[k][i][j] = 0d;
        } } }
        
        
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++)  {
	        fi[0][i][j] = 0d;
     	    for (int k = 1; k < 5; k++) {
   	            fi[k][i][j] = (fx[i][j]*cx[k] + fy[i][j]*cy[k])*rho[i][j]/3.0d;
            }
            for (int k = 5; k < 9; k++) {
                fi[k][i][j] = (fx[i][j]*cx[k] + fy[i][j]*cy[k])*rho[i][j]/12.0d;
            }
        } }
        
        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])/tau + fi[k][i][j];
        } } }
        
        
      	// 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];
        } } }
        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];
        } }
        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];
        } }
        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];
        } }
        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];
        } }
        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];
        } }
        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];
        } }
        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];
        } }
        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];
        } }

        //boundary condition
        for (int j = 0; j < ny; j++){
            f[1][   0][j] = f[3][   0][j];
            f[5][   0][j] = f[7][   0][j];
            f[8][   0][j] = f[6][   0][j];
            f[3][nx-1][j] = f[1][nx-1][j];
            f[7][nx-1][j] = f[5][nx-1][j];
            f[6][nx-1][j] = f[8][nx-1][j];
         }

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

        //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] = 0d; v[i][j] = 0d;
            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];
            }
            if(rho[i][j] != 0){
                u[i][j] = u[i][j]/rho[i][j];
                v[i][j] = v[i][j]/rho[i][j];
            }else{
                u[i][j] = 0d; v[i][j] = 0d;
            }
        } }

        // A circular ripple 
        if(tflag == 1){
            for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
                if(i >=0){ if(i <nx){ if(j >=0){ if(j <ny){ if((i-tx)*(i-tx)+(j-ty)*(j-ty) !=0){
                //solitary wave
                u[i][j] = 5.0d*(Math.exp(-(((i-tx)*(i-tx)+(j-ty)*(j-ty))/(4d*1d)))/(2d*Math.sqrt(Math.PI*1.0d))*(double)(i-tx)/(Math.sqrt((i-tx)*(i-tx)+(j-ty)*(j-ty))));
                v[i][j] = 5.0d*(Math.exp(-(((i-tx)*(i-tx)+(j-ty)*(j-ty))/(4d*1d)))/(2d*Math.sqrt(Math.PI*1.0d))*(double)(j-ty)/(Math.sqrt((i-tx)*(i-tx)+(j-ty)*(j-ty))));
                } } } } }
            } }
            tflag = 0;
        }
        tmp = 0f;
        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; }
        } }


         //S-function
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++){
            up[i][j] = 0.0d; vp[i][j] = 0.0d; phi[i][j] = 0.0d;
            for (int m = 0; m < np; m++) {
              tmp1 = rp - Math.sqrt((xp[m] - (double)(i))*(xp[m] - (double)(i)) + (yp[m] - (double)(j))*(yp[m] - (double)(j)));
              if(tmp1 > guzai/2.0d){
                tmp2 = 1.0d;
              }else if(tmp1 < - guzai/2.0d){
                tmp2 = 0.0d;
              }else{
                tmp2 = (Math.sin(Math.PI*tmp1/guzai) + 1.0d)/2.0d;
              }
              up[i][j] =  up[i][j] + tmp2*(upt[m] - omega[m]*((double)(j) - yp[m]));
              vp[i][j] =  vp[i][j] + tmp2*(vpt[m] + omega[m]*((double)(i) - xp[m]));
             phi[i][j] = phi[i][j] + tmp2;
            }
        } }

        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++){
   	    if(phi[i][j] != 0){
   	        up[i][j] = up[i][j]/phi[i][j];
   	        vp[i][j] = vp[i][j]/phi[i][j];
   	     }else{
   	        up[i][j] = 0.0d;
   	        vp[i][j] = 0.0d;
   	     }
        } }

        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++){
            fx[i][j] = phi[i][j]*(up[i][j] - u[i][j]);
            fy[i][j] = phi[i][j]*(vp[i][j] - v[i][j]);
        } }

        // particle-particle collisions
        for (int m = 0; m < np; m++) {
           fxc[m] = 0.0d; fyc[m] = 0.0d;
        }

        for (int m = 0; m < np; m++) { for (int n = 0; n < np; n++) {
            if( m != n){
              tmp1 = Math.sqrt((xp[m] - xp[n])*(xp[m] - xp[n]) + (yp[m] - yp[n])*(yp[m] - yp[n]));
              if(tmp1 < 2d*rp + zetap){
	        fxc[m] = fxc[m] + (xp[m] - xp[n])/Math.abs(xp[m] - xp[n])*(2d*rp - tmp1 + zetap)/zetap
                                                                         *(2d*rp - tmp1 + zetap)/zetap*0.5d;
	        fyc[m] = fyc[m] + (yp[m] - yp[n])/Math.abs(yp[m] - yp[n])*(2d*rp - tmp1 + zetap)/zetap
	                                                                 *(2d*rp - tmp1 + zetap)/zetap*0.5d;
              }
            }
        } }         
        
        // particle-wall collisions
        for (int m = 0; m < np; m++) {
          if (xp[m] < rp + zetaw) {
            fxc[m] = fxc[m] +xp[m]/Math.abs(xp[m])*(rp + zetaw - xp[m])/zetaw
                                                  *(rp + zetaw - xp[m])/zetaw*0.1d;
          }
          if (yp[m] < rp + zetaw) {
            fyc[m] = fyc[m] + yp[m]/Math.abs(yp[m])*(rp + zetaw - yp[m])/zetaw
                                                   *(rp + zetaw - yp[m])/zetaw*0.1d;
          }
          if (nx - xp[m] < rp + zetaw) {
            fxc[m] = fxc[m] +(xp[m] - nx)/Math.abs(xp[m] - nx)*(rp + zetaw - nx + xp[m])/zetaw
                                                              *(rp + zetaw - nx + xp[m])/zetaw*0.1d;
          }
          if (ny - yp[m] < rp + zetaw) {
            fyc[m] = fyc[m] + (yp[m] - ny)/Math.abs(yp[m] - ny)*(rp + zetaw - ny + yp[m])/zetaw
                                                               *(rp + zetaw - ny + yp[m])/zetaw*0.1d;
          }
        }
   
        // particle motion 
        for (int m = 0; m < np; m++) {
          fxp[m] = 0.0d; fyp[m] = 0.0d; torqp[m] = 0.0d;
          for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++){
            if((xp[m] - (double)(i))*(xp[m] - (double)(i)) + (yp[m] - (double)(j))*(yp[m] - (double)(j))
             <= (rp + guzai/2.0d)*(rp + guzai/2.0d) ){
   	        fxp[m] = fxp[m] + rho[i][j]*phi[i][j]*(u[i][j] - up[i][j]);
   	        fyp[m] = fyp[m] + rho[i][j]*phi[i][j]*(v[i][j] - vp[i][j]);
   	        torqp[m] =torqp[m] + ((double)(i) - xp[m])*fy[i][j] - ((double)(j) - yp[m])*fx[i][j];
            }
          }
        } }

        for (int m = 0; m < np; m++) {
          upt[m] = (1.0d + 1.0d/rhop)*up1[m] - 1.0d/rhop*up2[m]
                 + (fxp[m] + fxc[m])/mp + (1.0d - 1.0d/rhop)*gxp;
          vpt[m] = (1.0d + 1.0d/rhop)*vp1[m] - 1.0d/rhop*vp2[m]
                 + (fyp[m] + fyc[m])/mp + (1.0d - 1.0d/rhop)*gyp;
        omega[m] = (1.0d + 1.0d/rhop)*omega1[m] - 1.0d/rhop*omega2[m] - torqp[m]/iip;
           xp[m] =  xp[m] + (upt[m]   +    up1[m])*0.5d;
           yp[m] =  yp[m] + (vpt[m]   +    vp1[m])*0.5d;
        theta[m] = theta[m] + (omega[m] + omega1[m])*0.5d;

        up2[m] =    up1[m];    up1[m] =   upt[m];
        vp2[m] =    vp1[m];    vp1[m] =   vpt[m];
        omega2[m] = omega1[m]; omega1[m] = omega[m];
        }
 	    
        time += dt;

        //graphics
      	c.drawColor(Color.WHITE);

        //velocity
      	p.setStrokeWidth(2);
        for (int i = 0; i < nx; i++) { for (int j = 0; j < ny; j++) {
            u2 = Math.sqrt(u[i][j]*u[i][j] + v[i][j]*v[i][j]);
     	    p.setColor(Color.argb((int)(u2/u0*255),0,0,(int)(u2/u0*255)));
            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);
   	    p.setColor(Color.argb(255,0, 255, 255));
            c.drawLine((float)(x[i][j] + wid/2.0d),(float)(y[i][j] + wid/2.0d),
                       (float)(x[i][j] + wid/2.0d + u[i][j]/u0*wid),
                       (float)(y[i][j] + wid/2.0d + v[i][j]/u0*wid), p);
        } }

        //particle
        for (int m = 0; m < np; m++) {
          p.setStyle(Paint.Style.FILL);
          p.setColor(Color.argb(127,255*(np-m),127,127));
          c.drawCircle((float)(xp[m]*wid), (float)(yp[m]*wid), (float)(rp*wid), p);
          p.setColor(Color.WHITE);
          c.drawLine((float)((xp[m]-rp*Math.cos(theta[m]+Math.PI/2.0d))*wid),(float)((yp[m]-rp*Math.sin(theta[m]+Math.PI/2.0d))*wid),
                     (float)((xp[m]+rp*Math.cos(theta[m]+Math.PI/2.0d))*wid),(float)((yp[m]+rp*Math.sin(theta[m]+Math.PI/2.0d))*wid), p);
        }
        
        // Text
        t.setTextSize (12);
        t.setColor(Color.BLACK);
        String msg1 = "time: " + time;
        c.drawText(msg1, 10, 480, t);
        String msg2 = "Particle1 ( "+(float)xp[0] +" , " +(float)yp[0]+" )"+"Angle ( "+(float)theta[0]+" )";
        c.drawText(msg2, 10, 500, t);
        String msg3 = "Particle2 ( "+(float)xp[1] +" , " +(float)yp[1]+" )"+"Angle ( "+(float)theta[1]+" )";
        c.drawText(msg3, 10, 520, t);

        t.setTextSize (11);
        String msg4 = "Android Version of the Lattice Boltzmann Method combined with smoothed-profile method";
        c.drawText(msg4, 10, 540, t);
        String msg5 = "Copyright @2014- Takeshi Seta, Kento Tanaka, Masahiro Matsuzuru, Takayuki Yamagishi";
        c.drawText(msg5, 10, 560, t);
        String msg6 = "All Rights Reserved.";
        c.drawText(msg6, 10, 580, t);

      	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);
    }
}