· 

格子ボルツマン法による流体シミュレーション

(著)山たー

Processingで書いてみたもの。Processing.jsを使っているので、参考にしたFluid Dynamics Simulationに比べるととても遅い。プログラムの解説とアルゴリズムはいずれ書きますが、とりあえずプログラムだけ。

 

流体シミュレーション

プログラムの実行/停止


物理量の表示切替


障害物の表示切替


流線の表示切替

コードと解説

全体のコード

//計算のための変数・配列の設定
//表示に関する変数・配列
//x,yの次元数
int xdim=200;
int ydim=80;
//シミュレーション画面のサイズ
int DispWidth=1000;
int DispHeight=400;
int xedge,yedge;  //格子サイズ
int step=0; //計算のステップ数
int visc_x=837;//粘度のスクロールバー
int flow_x=850;//流速のスクロールバー
//計算シークエンス変数
int bar_seq=1;       //バリアの種類(0~7)
int disp_seq=1;      //画面表示(0~5)
int disp_flowVec=1;  //流線の表示(0 or 1)
int free_bar_seq=1;  //バリアの手書きモード(0 or 1)

//初期物理量の設定
float u=0.075;                 //初期流速
float usq=u*u;                    //初期流速の2乗
float viscosity=0.05;            //粘度
float omega=1/(3*viscosity+0.5); //単一時間緩和係数

//速度分布関数
float[][] f0=new float[xdim][ydim];
float[][] fN=new float[xdim][ydim];
float[][] fS=new float[xdim][ydim];
float[][] fE=new float[xdim][ydim];
float[][] fW=new float[xdim][ydim];
float[][] fNW=new float[xdim][ydim];
float[][] fNE=new float[xdim][ydim];
float[][] fSW=new float[xdim][ydim];
float[][] fSE=new float[xdim][ydim];

//分布関数より計算される物理量
float[][] density=new float[xdim][ydim]; //密度
float[][] xvel=new float[xdim][ydim];    //速度のx方向
float[][] yvel=new float[xdim][ydim];    //速度のy方向
float[][] speed2=new float[xdim][ydim];  //流速
float[][] rot=new float[xdim][ydim];     //回転(rotation)

//障害物設定のための変数・配列
//格子内の障壁の有無(有=1,無=0)
int[][] barrier=new int[xdim][ydim];
//円形(半径6)の障壁の位置座標の逆
int[] I_sircle_x={1,1,1,1,2,2,3,4};
int[] I_sircle_y={1,2,3,4,1,2,1,1};
int num_Isx=I_sircle_x.length;
//円形(半径20)の障壁の位置座標の逆
int[] I_big_sircle_x={1,1,1,1,1,1,1,1,1, 1, 1, 1, 1, 1, 1, 1,2,2,2,2,2,2,2,2,2, 2, 2, 2,
3,3,3,3,3,3,3,3,3, 3,4,4,4,4,4,4,4,4,4,5,5,5,5,5,5,5,6,6,6,6,6,6,7,7,7,7,7,8,8,8,8,9,9,9,9,10,10,10,
11,11,12,12,13,14,15,16};
int[] I_big_sircle_y={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,1,2,3,4,5,6,7,8,9,10,11,12,
1,2,3,4,5,6,7,8,9,10,1,2,3,4,5,6,7,8,9,1,2,3,4,5,6,7,1,2,3,4,5,6,1,2,3,4,5,1,2,3,4,1,2,3,4,1,2,3,
1,2,1,2,1,1,1,1};
int num_Ibsx=I_big_sircle_x.length;
//計算のショートカット
float four9ths = 4.0 / 9;
float one9th = 1.0 / 9;
float one36th = 1.0 / 36;

//メイン関数
void setup(){
  size(1000,700);
  colorMode(HSB);
  noStroke();
  frameRate(500);
  xedge=DispWidth/xdim;
  yedge=DispHeight/ydim;
  initBarrier(); //障害物の初期設定
  initFluid();   //流れの初期設定
}
void draw(){
  background(0);
  //計算
  Collision(); //衝突
  Streaming();  //並進
  BounceBack();  //境界
  
  //表示
  if(disp_seq==0){  
    fill(255);
    noStroke();
    rect(0,0,DispWidth,DispHeight);
  }else if(disp_seq==1){
    DispXvel();     //流速のx方向
  }else if(disp_seq==2){
    DispYvel();     //流速のy方向
  }else if(disp_seq==3){
    DispSpeed2();   //速さの2乗
  }else if(disp_seq==4){
    DispDensity();  //密度
  }else if(disp_seq==5){
    DispRotation();     //回転
  }
  if(disp_flowVec==1){
    DispFlowVector();  //流速ベクトル
  }
  DispBarrier(); //障害物の表示
  DispUI();      //UserInterfaceの表示
  step+=1;
  //println("FrameRate : "+int(frameRate) + ", Step : "+step);
}

//初期設定
//流れの初期設定
void initFluid(){
  step=0;
  for (int x=0; x<xdim; x++) {
    for (int y=0; y<ydim; y++) {
      if (barrier[x][y]==1) {
        zeroLattice(x,y);
      }else{
        //usq=u^2
        f0[x][y]=four9ths*(1-1.5*usq);
        fN[x][y]=one9th*(1-1.5*usq);
        fS[x][y]=one9th*(1-1.5*usq);
        fE[x][y]=one9th*(1+3*u+3*usq);
        fW[x][y]=one9th*(1-3*u+3*usq);
        fNE[x][y]=one36th*(1+3*u+3*usq);
        fSE[x][y]=one36th*(1+3*u+3*usq);
        fNW[x][y]=one36th*(1-3*u+3*usq);
        fSW[x][y]=one36th*(1-3*u+3*usq);
        density[x][y]=1;
        xvel[x][y]=u;
        yvel[x][y]=0;
        speed2[x][y]=usq;
      }
    }
  }
}
//流れがない場
void zeroLattice(int x,int y){
  f0[x][y]=0;
  fE[x][y]=0;
  fW[x][y]=0;
  fN[x][y]=0;
  fS[x][y]=0;
  fNE[x][y]=0;
  fNW[x][y]=0;
  fSE[x][y]=0;
  fSW[x][y]=0;
  xvel[x][y]=0;
  yvel[x][y]=0;
  speed2[x][y]=0;
}
//障害物の設定
void initBarrier(){
  if(bar_seq==0){
     clearBarrier();
  }else if(bar_seq==1){
    for(int x=50;x<70;x++){
      for(int y=ydim-20;y<ydim-1;y++){
        barrier[x][y]=1;
      }
    }
  }else if(bar_seq==2){
    for(int x=51;x<91;x++){
      for(int y=20;y<61;y++){
        barrier[x][y]=1;
      }
    }
    for(int x=51;x<88;x++){
      for(int y=23;y<58;y++){
        barrier[x][y]=0;
      }
    }
  }else if(bar_seq==3){
    for(int x=51;x<63;x++){
      for(int y=35;y<47;y++){
        barrier[x][y]=1;
      }
    }
    for(int i=0;i<num_Isx;i++){
      barrier[50+I_sircle_x[i]][34+I_sircle_y[i]]=0;
      barrier[50+(13-I_sircle_x[i])][34+I_sircle_y[i]]=0;
      barrier[50+I_sircle_x[i]][34+(13-I_sircle_y[i])]=0;
      barrier[50+(13-I_sircle_x[i])][34+(13-I_sircle_y[i])]=0;
    }
  }else if(bar_seq==4){
    for(int x=31;x<71;x++){
      for(int y=22;y<62;y++){
        barrier[x][y]=1;
      }
    }
    for(int i=0;i<num_Ibsx;i++){
      barrier[30+I_big_sircle_x[i]][21+I_big_sircle_y[i]]=0;
      barrier[30+(41-I_big_sircle_x[i])][21+I_big_sircle_y[i]]=0;
      barrier[30+I_big_sircle_x[i]][21+(41-I_big_sircle_y[i])]=0;
      barrier[30+(41-I_big_sircle_x[i])][21+(41-I_big_sircle_y[i])]=0;
    }
  }else if(bar_seq==5){
    for(int x=40;x<43;x++){
      for(int y=20;y<61;y++){
        barrier[x][y]=1;
      }
    }
  }else if(bar_seq==6){
    for(int x=51;x<101;x++){
      for(int y=39;y<41;y++){
        barrier[x][y]=1;
      }
    }
  }else if(bar_seq==7){
    for(int x=31;x<51;x++){
        barrier[100-x][x]=1;
      }
  }
}
//全体の障害物を消去
void clearBarrier(){
  for(int x=0;x<xdim;x++){
    for(int y=0;y<ydim;y++){
      barrier[x][y]=0;
    }
  }
}

//LBM計算関数
//各粒子の衝突
void Collision() {
  float rho,Irho,one9thRho,one36thRho,vx,vy,vxx,vyy,vx3,vy3,vxy,vsq,vsq15;
  for (int x=0; x<xdim; x++) {
    for (int y=0; y<ydim; y++) {
      //障害物がない場合
      if (barrier[x][y]==0) {
        //格子内の全粒子数=密度ρ
        rho=f0[x][y]+fN[x][y]+fS[x][y]+fE[x][y]+fW[x][y]+fNW[x][y]+fNE[x][y]+fSW[x][y]+fSE[x][y];
        //密度(macroscopic density)(単位空間内なので個数に等しい)
        density[x][y]=rho;
        one9thRho=one9th*rho; //ρ/9
        one36thRho=one36th*rho; //ρ/36
        Irho=1/rho; ///inverse rho :ρの逆数
        if(rho>0){
          vx=(fE[x][y]+fNE[x][y]+fSE[x][y]-fW[x][y]-fNW[x][y]-fSW[x][y])*Irho;
          vy=(fN[x][y]+fNE[x][y]+fNW[x][y]-fS[x][y]-fSE[x][y]-fSW[x][y])*Irho;
        }else{
          vx = 0;
          vy = 0;
        }
        xvel[x][y]=vx; //x方向の速度
        yvel[x][y]=vy; //y方向の速度
        
        vx3=3*vx;
        vy3=3*vy;
        vxx=vx*vx;
        vyy=vy*vy;
        vxy=2*vx*vy;
        vsq=vxx+vyy;
        speed2[x][y]=vsq;
        vsq15=1.5*vsq;
        f0[x][y] +=omega*(four9ths*rho *(1                       -vsq15)-f0[x][y]);
        fE[x][y] +=omega*(   one9thRho *(1+vx3     +4.5*vxx      -vsq15)-fE[x][y]);
        fW[x][y] +=omega*(   one9thRho *(1-vx3     +4.5*vxx      -vsq15)-fW[x][y]);
        fN[x][y] +=omega*(   one9thRho *(1+vy3     +4.5*vyy      -vsq15)-fN[x][y]);
        fS[x][y] +=omega*(   one9thRho *(1-vy3     +4.5*vyy      -vsq15)-fS[x][y]);
        fNE[x][y]+=omega*(  one36thRho *(1+vx3+ vy3+4.5*(vsq+vxy)-vsq15)-fNE[x][y]);
        fNW[x][y]+=omega*(  one36thRho *(1-vx3+ vy3+4.5*(vsq-vxy)-vsq15)-fNW[x][y]);
        fSE[x][y]+=omega*(  one36thRho *(1+vx3- vy3+4.5*(vsq-vxy)-vsq15)-fSE[x][y]);
        fSW[x][y]+=omega*(  one36thRho *(1-vx3- vy3+4.5*(vsq+vxy)-vsq15)-fSW[x][y]);
      }
    }
  }
}
//粒子の移動
void Streaming(){
  for (int x=0; x<xdim-1; x++) { 
    for (int y=ydim-1; y>0; y--) {
      fN[x][y] =fN[x][y-1];   
      fNW[x][y]=fNW[x+1][y-1]; 
    }
  }
  for (int x=xdim-1; x>0; x--) {   
    for (int y=ydim-1; y>0; y--) {
      fE[x][y] =fE[x-1][y];   
      fNE[x][y]=fNE[x-1][y-1]; 
    }
  }
  for (int x=xdim-1; x>0; x--) {   
    for (int y=0; y<ydim-1; y++) {
      fS[x][y] =fS[x][y+1];  
      fSE[x][y]=fSE[x-1][y+1]; 
    }
  }
  for (int x=0; x<xdim-1; x++) {  
    for (int y=0; y<ydim-1; y++) {
      fW[x][y] =fW[x+1][y];   
      fSW[x][y]=fSW[x+1][y+1];  
    }
  }
  for (int y=0; y<ydim-1; y++) {
    fS[0][y]=fS[0][y+1];
  }
  for (int y=ydim-1; y>0; y--) {
    fN[xdim-1][y]=fN[xdim-1][y-1];
  }
  
  //左端から新規に流入させる
  for (int y=0; y<ydim; y++) {
    if (barrier[0][y]==0) {
      fE[0][y] =one9th *(1+3*u+3*usq);
      fNE[0][y]=one36th*(1+3*u+3*usq);
      fSE[0][y]=fNE[0][y];
      //fSE[0][y]=one36th*(1+3*u+3*usq);
    }
  }
  //上下の境界は初期の流れと同じにしておく=計算外にするということ
  for (int x=0; x<xdim; x++) {
    f0[x][0] =four9ths*(1-1.5*usq);
    fE[x][0] =  one9th*(1+3*u+3*usq);
    fW[x][0] =  one9th*(1-3*u+3*usq);
    fN[x][0] =  one9th*(1-1.5*usq);
    fS[x][0] =  one9th*(1-1.5*usq);
    fNE[x][0]= one36th*(1+3*u+3*usq);
    fSE[x][0]= one36th*(1+3*u+3*usq);
    fNW[x][0]= one36th*(1-3*u+3*usq);
    fSW[x][0]= one36th*(1-3*u+3*usq);
    f0[x][ydim-1]=four9ths*(1-1.5*usq);
    fE[x][ydim-1]=  one9th*(1+3*u+3*usq);
    fW[x][ydim-1]=  one9th*(1-3*u+3*usq);
    fN[x][ydim-1]=  one9th*(1-1.5*usq);
    fS[x][ydim-1]=  one9th*(1-1.5*usq);
    fNE[x][ydim-1]= one36th*(1+3*u+3*usq);
    fSE[x][ydim-1]= one36th*(1+3*u+3*usq);
    fNW[x][ydim-1]= one36th*(1-3*u+3*usq);
    fSW[x][ydim-1]= one36th*(1-3*u+3*usq);
  }
  //流出端は1つ左隣の分布と同じとする(これがないと不安定化する)
  for(int y=0; y<ydim; y++){
    if (barrier[xdim-1][y]==0){
      fW[xdim-1][y]=fW[xdim-2][y];
      fNW[xdim-1][y]=fNW[xdim-2][y];
      fSW[xdim-1][y]=fSW[xdim-2][y];
    }
  }
}
//跳ね返り
void BounceBack(){
  //Bounce-back条件, 単に180°回転させる
  for (int x=0; x<xdim; x++) {
    for (int y=0; y<ydim; y++) {
      if (barrier[x][y]==1) {
        if (fN[x][y]>0){fS[x][y-1]+=fN[x][y]; fN[x][y]=0; }
        if (fS[x][y]>0){fN[x][y+1]+=fS[x][y]; fS[x][y]=0; }
        if (fE[x][y]>0){fW[x-1][y]+=fE[x][y]; fE[x][y]=0; }
        if (fW[x][y]>0){fE[x+1][y]+=fW[x][y]; fW[x][y]=0; }
        if (fNW[x][y]>0){fSE[x+1][y-1]+=fNW[x][y]; fNW[x][y]=0; }
        if (fNE[x][y]>0){fSW[x-1][y-1]+=fNE[x][y]; fNE[x][y]=0; }
        if (fSW[x][y]>0){fNE[x+1][y+1]+=fSW[x][y]; fSW[x][y]=0; }
        if (fSE[x][y]>0){fNW[x-1][y+1]+=fSE[x][y]; fSE[x][y]=0; }
      }
    }
  }
}

//サブ計算・結果表示関数
//回転(rotation)の計算
void ComputeRotation(){
  for (int x=1; x<xdim-1; x++) {
    for (int y=1; y<ydim-1; y++) {
      rot[x][y]=(yvel[x+1][y]-yvel[x-1][y])-(xvel[x][y+1]-xvel[x][y-1]);
    }
  }
  //端点は気を付ける
  for (int y=1; y<ydim-1; y++) {
    rot[0][y]=2*(yvel[1][y]-yvel[0][y])-(xvel[0][y+1]-xvel[0][y-1]);
    rot[xdim-1][y]=2*(yvel[xdim-1][y]-yvel[xdim-2][y])-(xvel[xdim-1][y+1]-xvel[xdim-1][y-1]);
  }
}
//rotの計算と表示
void DispRotation(){
  ComputeRotation();
  noStroke();
  for(int x=1;x<xdim-1;x++){
    for(int y=1;y<ydim-1;y++){
      int col=constrain(int(rot[x][y]*4000),-255,255);
      if(col>0){
        fill(0,col,255);
      }else if(col<-0){
        fill(160,-col,255);
      }else if(col==0){
        fill(255);
      }
      rect(xedge*x,yedge*y,xedge,yedge);
    }
  }
}
//x方向の流速の表示
void DispXvel(){
  noStroke();
  for(int x=1;x<xdim-1;x++){
    for(int y=1;y<ydim-1;y++){
      int col=constrain(int((xvel[x][y]*20)*100),0,220);
      fill(col,255,255);
      rect(xedge*x,yedge*y,xedge,yedge);
    }
  }
}
//y方向の流速の表示
void DispYvel(){
  noStroke();
  for(int x=1;x<xdim-1;x++){
    for(int y=1;y<ydim-1;y++){
      int col=constrain(int((yvel[x][y]*20)*100),0,220);
      fill(col,255,255);
      rect(xedge*x,yedge*y,xedge,yedge);
    }
  }
}
//流速の表示
void DispSpeed2(){
  noStroke();
  for(int x=1;x<xdim-1;x++){
    for(int y=1;y<ydim-1;y++){
      int col=constrain(int((speed2[x][y]*20)*100),0,220);
      fill(col,255,255);
      rect(xedge*x,yedge*y,xedge,yedge);
    }
  }
}
//密度の表示
void DispDensity(){
  noStroke();
  for(int x=1;x<xdim-1;x++){
    for(int y=1;y<ydim-1;y++){
      int col=int((0.5+density[x][y]*100*0.2)*10);
      fill(col,255,255);
      rect(xedge*x,yedge*y,xedge,yedge);
    }
  }
}
//流速ベクトルの表示
void DispFlowVector(){
  stroke(0);
  for(int x=1;x<xdim-1;x++){
    if(x%5==0){
      float cx=xedge*x+xedge/2;
      for(int y=1;y<ydim-1;y++){
        if(x%5==0 && y%5==0){
          float cy=yedge*y+yedge/2;
          line(cx-xvel[x][y]*50,cy-yvel[x][y]*50,cx+xvel[x][y]*50,cy+yvel[x][y]*50);
        }
      }
    }
  }
}
//障害物の表示
void DispBarrier(){
 for(int x=1;x<xdim-1;x++){
  for(int y=1;y<ydim-1;y++){
   if(barrier[x][y]==1){
      fill(0);
      rect(xedge*x,yedge*y,xedge,yedge);
     }
    }
  }
}

//User Interfaceの表示
void DispUI(){
  fill(255);
  textSize(20);
  textAlign(CENTER);
  text("Fluid Dynamics Simulation with Lattice Boltzmann Method", width/2, DispHeight+25);
  
  //Display Style選択の表示
  textSize(25);
  textAlign(LEFT);
  text("Display Style",70, DispHeight+60);
  textSize(20);
  text("X Velocity",120, DispHeight+90);
  text("Y Velocity",120, DispHeight+120);
  text("Speed",120, DispHeight+150);
  text("Density",120, DispHeight+180);
  text("Rotation",120, DispHeight+210);
  text("None",120, DispHeight+240);
  text("FlowVector",120, DispHeight+280);
  
  //ボタンの作成
  stroke(0);
  rect(85,DispHeight+70,22,22);
  rect(85,DispHeight+100,22,22);
  rect(85,DispHeight+130,22,22);
  rect(85,DispHeight+160,22,22);
  rect(85,DispHeight+190,22,22);
  rect(85,DispHeight+220,22,22);
  rect(85,DispHeight+260,22,22);
  
  //選択中のモードの表示
  fill(130,255,255);
  if(disp_seq==0){  
    rect(87,DispHeight+222,18,18);
  }else if(disp_seq==1){
    rect(87,DispHeight+72,18,18);
  }else if(disp_seq==2){
    rect(87,DispHeight+102,18,18);
  }else if(disp_seq==3){
    rect(87,DispHeight+132,18,18);
  }else if(disp_seq==4){
    rect(87,DispHeight+162,18,18);
  }else if(disp_seq==5){
    rect(87,DispHeight+192,18,18);
  }
  if(disp_flowVec==1){
    rect(87,DispHeight+262,18,18);
  }
  
  //障害物選択の表示
  fill(255);
  textSize(25);
  textAlign(LEFT);
  text("Barrier",320, DispHeight+60);
  textSize(20);
  text("Box1",370, DispHeight+90);
  text("Box2",370, DispHeight+120);
  text("Circle",370, DispHeight+150);
  text("Big Circle",370, DispHeight+180);
  text("Board1",370, DispHeight+210);
  text("Board2",370, DispHeight+240);
  text("Board3",370, DispHeight+270);
  text("None",520, DispHeight+90);
  text("Freely Writing",520, DispHeight+130);
  stroke(0);
  rect(335,DispHeight+70,22,22);
  rect(335,DispHeight+100,22,22);
  rect(335,DispHeight+130,22,22);
  rect(335,DispHeight+160,22,22);
  rect(335,DispHeight+190,22,22);
  rect(335,DispHeight+220,22,22);
  rect(335,DispHeight+250,22,22);
  rect(485,DispHeight+70,22,22);
  rect(485,DispHeight+110,22,22);
  
  //選択中のモードの表示
  fill(0,255,255);
  if(bar_seq==0){  
    rect(487,DispHeight+72,18,18);
  }else if(bar_seq==1){
    rect(337,DispHeight+72,18,18);
  }else if(bar_seq==2){
    rect(337,DispHeight+102,18,18);
  }else if(bar_seq==3){
    rect(337,DispHeight+132,18,18);
  }else if(bar_seq==4){
    rect(337,DispHeight+162,18,18);
  }else if(bar_seq==5){
    rect(337,DispHeight+192,18,18);
  }else if(bar_seq==6){
    rect(337,DispHeight+222,18,18);
  }else if(bar_seq==7){
    rect(337,DispHeight+252,18,18);
  }
  if(free_bar_seq==1){  
    rect(487,DispHeight+112,18,18);
  }
  
  //粘度の変更
  fill(255);
  textSize(25);
  textAlign(LEFT);
  text("Viscosity",720, DispHeight+60);
  text(viscosity,870, DispHeight+60);
  stroke(255);
  line(750,DispHeight+100,970,DispHeight+100);
  line(750,DispHeight+90,750,DispHeight+110);
  line(970,DispHeight+90,970,DispHeight+110);
  textAlign(CENTER);
  textSize(15);
  text("0.005",750,DispHeight+130);
  text("1",970,DispHeight+130);
  fill(100,255,255,200);
  rect(visc_x, DispHeight+85, 20, 30);
  
  //流速の変更
  fill(255);
  textSize(25);
  textAlign(LEFT);
  text("Flow Speed",720, DispHeight+180);
  text(u,870, DispHeight+180);
  stroke(255);
  line(750,DispHeight+220,970,DispHeight+220);
  line(750,DispHeight+210,750,DispHeight+230);
  line(970,DispHeight+210,970,DispHeight+230);
  textAlign(CENTER);
  textSize(15);
  text("0.05",750,DispHeight+250);
  text("0.10",970,DispHeight+250);
  fill(150,255,255,200);
  rect(flow_x, DispHeight+205, 20, 30);
}

//キーボード・マウスとのインタラクション
//ドラッグされたとき
void mouseDragged(){
  if(free_bar_seq==1){
    if(mouseY<DispHeight){
      int x=constrain(mouseX,10,DispWidth-20)/xedge;
      int y=constrain(mouseY,10,DispHeight-10)/yedge;
      barrier[x][y]=1;
    }
  }
  if(mouseX>=750 && mouseX<=950 && mouseY>=DispHeight+90 && mouseY<=DispHeight+110){
    visc_x=mouseX;
    viscosity=exp((visc_x-750)*0.0265-5.3);
    omega=1/(3*viscosity+0.5);
  }
  if(mouseX>=750 && mouseX<=950 && mouseY>=DispHeight+210 && mouseY<=DispHeight+230){
    flow_x=mouseX;
    u=(flow_x-750)*0.00025+0.05;
    usq=u*u;
  }
}
//クリックされたとき
void mouseClicked(){
  if(free_bar_seq==1){
    if(mouseY<DispHeight){
      int x=constrain(mouseX,10,DispWidth-20)/xedge;
      int y=constrain(mouseY,10,DispHeight-10)/yedge;
      barrier[x][y]=1;
    }
  }
  //表示モードの変更
  if(mouseX>85 && mouseX<107 && mouseY>DispHeight+70 && mouseY<DispHeight+92){
    disp_seq=1;
  }else if(mouseX>85 && mouseX<107 && mouseY>DispHeight+100 && mouseY<DispHeight+122){
    disp_seq=2;
  }else if(mouseX>85 && mouseX<107 && mouseY>DispHeight+130 && mouseY<DispHeight+152){
    disp_seq=3;
  }else if(mouseX>85 && mouseX<107 && mouseY>DispHeight+160 && mouseY<DispHeight+182){
    disp_seq=4;
  }else if(mouseX>85 && mouseX<107 && mouseY>DispHeight+190 && mouseY<DispHeight+212){
    disp_seq=5;
  }else if(mouseX>85 && mouseX<107 && mouseY>DispHeight+220 && mouseY<DispHeight+242){
    disp_seq=0;
  }
  if(mouseX>85 && mouseX<107 && mouseY>DispHeight+260 && mouseY<DispHeight+282){
    if(disp_flowVec==1){
      disp_flowVec=0;
    }else{
      disp_flowVec=1;
    }
  }
  
  //障害物の表示
  if(mouseX>335 && mouseX<357 && mouseY>DispHeight+70 && mouseY<DispHeight+92){
    bar_seq=1;
    clearBarrier();
    initBarrier(); 
    initFluid();
  }else if(mouseX>335 && mouseX<357 && mouseY>DispHeight+100 && mouseY<DispHeight+122){
    bar_seq=2;
    clearBarrier();
    initBarrier(); 
    initFluid();
  }else if(mouseX>335 && mouseX<357 && mouseY>DispHeight+130 && mouseY<DispHeight+152){
    bar_seq=3;
    clearBarrier();
    initBarrier(); 
    initFluid();
  }else if(mouseX>335 && mouseX<357 && mouseY>DispHeight+160 && mouseY<DispHeight+182){
    bar_seq=4;
    clearBarrier();
    initBarrier(); 
    initFluid();
  }else if(mouseX>335 && mouseX<357 && mouseY>DispHeight+190 && mouseY<DispHeight+212){
    bar_seq=5;
    clearBarrier();
    initBarrier(); 
    initFluid();
  }else if(mouseX>335 && mouseX<357 && mouseY>DispHeight+220 && mouseY<DispHeight+242){
    bar_seq=6;
    clearBarrier();
    initBarrier(); 
    initFluid();
  }else if(mouseX>335 && mouseX<357 && mouseY>DispHeight+250 && mouseY<DispHeight+272){
    bar_seq=7;
    clearBarrier();
    initBarrier(); 
    initFluid();
  }else if(mouseX>485 && mouseX<507 && mouseY>DispHeight+70 && mouseY<DispHeight+92){
    bar_seq=0;
    clearBarrier();
    initBarrier(); 
    initFluid();
  }
  if(mouseX>485 && mouseX<507 && mouseY>DispHeight+110 && mouseY<DispHeight+130){
    if(free_bar_seq==1){
      free_bar_seq=0;
    }else{
      free_bar_seq=1;
    }
  }
}
//cを押すと初期状態に戻す
void keyPressed() {
  if ( key == 'c' ) {
    initBarrier(); 
    initFluid();
  }
}
}