とりあえず貼り付けときますね。
「時間感覚」のところと「粒子が近づきすぎてないか」のところが追加した処理。
import java.awt.*; import java.awt.image.BufferedImage; import javax.swing.*; public class Particle { enum Type{ WATER(0, new Color(0, 0, 255, 96)), OUTWALL(1, Color.BLACK), INWALL(2, Color.GREEN); int index; Color color; Type(int index, Color c){ color = c; this.index = index; } } static class Elem{ double x; double y; double vx; double vy; double p; Type type; boolean surface = false; public Elem(double x, double y, Type type) { this.x = x; this.y = y; this.type = type; } Color getColor(){ return surface ? new Color(255, 255, 255, 96) : type.color; } } static Elem[] elements = new Elem[1500]; static int elementCount; static JLabel label; static Image img = new BufferedImage(500, 400, BufferedImage.TYPE_INT_RGB); static final double parcSize = 8e-3;//1.5e-2;//粒子サイズ static final double minlen = 0.5;//粒子の最低距離(粒子サイズに対する) static final double re = 2.1;//近傍粒子の範囲 static final double relap = 4.0;//ラプラシアン用の近傍粒子の範囲 static double deltat = 0.0001;//時間の単位(粒子の最高速度によって調整される) static final double grav = 9.8;//重力加速度(m/s^2) static final double dim = 2;//次元数 static final double nyu = 1.0038E-6;//動粘性係数(水 1.0038E-6 m^2/s) static final double rho = 1.0e3;//粒子の比重(kg/m^3) static final double deltatmax = 1.0e-4;//時間の単位の最大値(本来は1e-3) static double n0;//基準粒子密度 static double lambda;//ラプラシアン用の係数 static final double beta = 0.97;//自由表面条件 public static void main(String[] args){ //左右の外壁 double maxy = 0; double oldparc = 1.5e-2; double ysize = 0.5 * parcSize / oldparc; double xsize = 0.8 * parcSize / oldparc; for(double y = 0; y < ysize; y += parcSize){ elements[elementCount++] = new Elem(0, y, Type.INWALL); elements[elementCount++] = new Elem(-parcSize, y, Type.OUTWALL); elements[elementCount++] = new Elem(-parcSize * 2, y, Type.OUTWALL); elements[elementCount++] = new Elem(xsize, y, Type.INWALL); elements[elementCount++] = new Elem(xsize+parcSize, y, Type.OUTWALL); elements[elementCount++] = new Elem(xsize+parcSize * 2, y, Type.OUTWALL); maxy = y; } //下の壁 for(double x = -parcSize * 2; x < xsize + parcSize * 2; x += parcSize){ elements[elementCount++] = new Elem(x, maxy + parcSize, (x < 0 || x > xsize) ? Type.OUTWALL : Type.INWALL); elements[elementCount++] = new Elem(x, maxy + parcSize * 2, Type.OUTWALL); elements[elementCount++] = new Elem(x, maxy + parcSize * 3, Type.OUTWALL); } //水柱 int sample = 0; double waterx = .3 * parcSize / oldparc; double watercxmin = .2 * parcSize / oldparc; double watercxmax = .25 * parcSize / oldparc; double watercymin = .1 * parcSize / oldparc; double watercymax = .15 * parcSize / oldparc; for(double x = parcSize; x < waterx; x+= parcSize){ for(double y = parcSize * 3; y < ysize; y += parcSize){//.5まで if(y > watercxmin && y < watercxmax && x > watercymin && x < watercymax){ //内部になる点を一点記憶しておく sample = elementCount; } elements[elementCount++] = new Elem(x, y, Type.WATER); } } //基本の粒子密度n0とラプラシアン用の係数rambda n0 = 0; double a = 0; double b = 0; for(int i = 0; i < elementCount; ++i){ if(i == sample) continue; double d = Math.sqrt( (elements[sample].x - elements[i].x) * (elements[sample].x - elements[i].x) + (elements[sample].y - elements[i].y) * (elements[sample].y - elements[i].y)); double w = (d < re * parcSize) ? parcSize * re / d - 1 : 0; n0 += w; a += (elements[sample].x - elements[i].x) * (elements[sample].x - elements[i].x) + (elements[sample].y - elements[i].y) * (elements[sample].y - elements[i].y) * w; } lambda = a / n0; //ウィンドウの表示 JFrame f = new JFrame("粒子法"); f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE); f.setSize(500, 450); label = new JLabel(new ImageIcon(img)); draw(); f.add(label); f.setVisible(true); //処理スレッド new Thread(){ @Override public void run() { try{ for(;;){ calc(); draw(); } }catch(InterruptedException e){} } }.start(); } static void calc() throws InterruptedException{ final double[][] w = new double[elementCount][elementCount]; //粒子密度を計算しておく for(int i = 0; i < elementCount; ++i){ for(int j = 0; j < i; ++j){ double d = Math.sqrt( (elements[i].x - elements[j].x) * (elements[i].x - elements[j].x) + (elements[i].y - elements[j].y) * (elements[i].y - elements[j].y)); double wd = (d < relap * parcSize) ? parcSize * relap / d - 1 : 0; w[i][j] = wd; w[j][i] = wd; } } //時間間隔 double maxvelo = 0; double velolim = 0.2 * parcSize / (deltatmax * 0.02); double svelolim = velolim * velolim; for(int i = 0; i < elementCount; ++i){ if(elements[i].type != Type.WATER) continue; double v = elements[i].vx * elements[i].vx + elements[i].vy * elements[i].vy; if(v > velolim){ //速すぎる場合、遅くしちゃう double sv = Math.sqrt(v); elements[i].vx *= svelolim / sv; elements[i].vy *= svelolim / sv; v = velolim; } if(v > maxvelo) maxvelo = v; } maxvelo = Math.sqrt(maxvelo); double dtlim = deltat * 1.1; if(maxvelo == 0){ deltat = dtlim; }else{ deltat = 0.2 * parcSize / maxvelo; } if(deltat > dtlim) deltat = dtlim; if(deltat > deltatmax) deltat = deltatmax; //仮の速度・位置を計算 final double oldx[] = new double[elementCount]; final double oldy[] = new double[elementCount]; for(int i = 0; i < elementCount; ++i){ oldx[i] = elements[i].vx; oldy[i] = elements[i].vy; } for(int i = 0; i < elementCount; ++i){ if(elements[i].type != Type.WATER){ continue; } //速度のラプラシアンを計算 double laplasux = 0; double laplasuy = 0; for(int j = 0; j < elementCount; ++j){ if(j == i) continue; if(w[i][j] == 0) continue; laplasux += (oldx[j] - elements[i].vx) * w[i][j]; laplasuy += (oldy[j] - elements[i].vy) * w[i][j]; } laplasux *= 2 * dim / n0 / lambda; laplasuy *= 2 * dim / n0 / lambda; //仮の位置・速度を計算 // elements[i].vx += + deltat * (nyu * laplasux); elements[i].vy += deltat * (nyu * laplasuy + grav); elements[i].x += elements[i].vx * deltat; elements[i].y += elements[i].vy * deltat; } //粒子が近づきすぎてないか double r2lim = parcSize * parcSize * minlen * minlen; for(int i = 0; i < elementCount; ++i){ for(int j = 0; j < i; ++j){ if(w[i][j] == 0) continue; double x = elements[j].x - elements[i].x; double y = elements[j].y - elements[i].y; double rr = x * x + y * y; if(rr < r2lim){ double r = Math.sqrt(rr); //double m2 = rho + rho; double vgx = (elements[i].vx + elements[j].vx) / 2;//本来は速度を粘度で按分 double vgy = (elements[i].vy + elements[j].vy) / 2; double vrx = rho * (elements[i].vx - vgx); double vry = rho * (elements[i].vy - vgy); double vabs = (vrx * x + vry * y) / r; if(vabs < 0) continue; final double colrat = 0.15; double vrat = 1 + colrat; double vmx = vrat * vabs * x / r; double vmy = vrat * vabs * y / r; if(elements[i].type == Type.WATER){ elements[i].vx -= vmx / rho; elements[i].vy -= vmy / rho; elements[i].x -= deltat * vmx / rho; elements[i].y -= deltat * vmy / rho; } if(elements[j].type == Type.WATER){ elements[j].vx += vmx / rho; elements[j].vy += vmy / rho; elements[j].x += deltat * vmx / rho; elements[j].y += deltat * vmy / rho; } } } } //仮の粒子密度を計算 final double[][] wasta = new double[elementCount][elementCount]; final double[][] rsqu = new double[elementCount][elementCount]; final int[] neighborcount = new int[elementCount]; final int[][] neighbor = new int[elementCount][elementCount]; for(int i = 0; i < elementCount; ++i){ for(int j = 0; j < i; ++j){ rsqu[i][j] = rsqu[j][i] = (elements[i].x - elements[j].x) * (elements[i].x - elements[j].x) + (elements[i].y - elements[j].y) * (elements[i].y - elements[j].y); double d = Math.sqrt(rsqu[i][j]); double rp = re * parcSize; if(d < rp){ double wd = rp / d - 1; wasta[i][j] = wd; wasta[j][i] = wd; neighbor[i][neighborcount[i]++] = j; neighbor[j][neighborcount[j]++] = i; } } } //粒子密度から自由表面になるところを検出 final double[] nasta = new double[elementCount]; for(int i = 0; i < elementCount; ++i){ for(int j = 0; j < neighborcount[i]; ++j){ nasta[i] += wasta[i][neighbor[i][j]]; } if(elements[i].type == Type.WATER){ elements[i].surface = nasta[i] < n0 * beta; } } //よくわからんけど。 final int[] check = new int[elementCount]; for(int i = 0; i < elementCount; ++i){ if(elements[i].type != Type.WATER){ check[i] = 2; }else if(elements[i].surface){ check[i] = 1; }else{ check[i] = 0; } } for(boolean cont = true; !cont;){ cont = false; for(int i = 0; i < elementCount; ++i){ if(check[i] == 1){ for(int j = 0; j < neighborcount[i]; ++j){ if(check[neighbor[i][j]] == 0) check[neighbor[i][j]] = 1; } check[i] = 2; cont = true; } } } //連立方程式の係数行列を作成 final double[] rk = new double[elementCount]; final double[][] coef = new double[elementCount][elementCount]; final double[] val = new double[elementCount]; //exec = Executors.newFixedThreadPool(threadLimit); for(int i = 0; i < elementCount; ++i){ if(elements[i].surface){ val[i] = 0; }else{ val[i] = -rho / deltat / deltat * (nasta[i] - n0) * lambda / 2 / dim; } for(int j = 0; j < neighborcount[i]; ++j){ coef[i][j + 1] = wasta[i][neighbor[i][j]]; coef[i][0] -= wasta[i][neighbor[i][j]]; } coef[i][0] -= 1.0e-7 / deltat / deltat * rho * n0;//よくわからないけどソースから if(check[i] == 0) coef[i][0] *= 2; } //連立方程式を解く(とりあえずガウスザイデル法。CG法がいいらしい) for(int n = 0; n < 1000; ++n){ double div = 0; for(int i = 0; i < elementCount; ++i){ double sum = 0; for(int j = 0; j < neighborcount[i]; ++j){ sum += coef[i][j + 1] * rk[neighbor[i][j]]; } double old = rk[i]; rk[i] = (val[i] - sum) / coef[i][0]; div += Math.abs(rk[i] - old); } //収束判定 if(div < 0.0001) { break; } } for(int i = 0; i < elementCount; ++i){ //圧力が負の場合は0に if(rk[i] < 0) rk[i] = 0; } //速度・場所を再計算 for(int i = 0; i < elementCount; ++i){ oldx[i] = elements[i].x; oldy[i] = elements[i].y; } for(int i = 0; i < elementCount; ++i){ final int i = idx; if(elements[i].type != Type.WATER) continue;//場所がかわるのは水だけ double udashx = 0; double udashy = 0; //圧力の最小値を得る double hatp = rk[i]; for(int j = 0; j < neighborcount[i]; ++j){ if(hatp > rk[neighbor[i][j]]) hatp = rk[neighbor[i][j]]; } //速度の修正量を計算 for(int j = 0; j < neighborcount[i]; ++j){ int n = neighbor[i][j]; udashx += (rk[n] - hatp) / rsqu[i][n] * (oldx[n] - elements[i].x) * wasta[i][n]; udashy += (rk[n] - hatp) / rsqu[i][n] * (oldy[n] - elements[i].y) * wasta[i][n]; } udashx = -udashx * deltat / rho * dim / n0; udashy = -udashy * deltat / rho * dim / n0; //速度の再計算 elements[i].vx += udashx; elements[i].vy += udashy; //場所を再計算 elements[i].x += deltat * udashx; elements[i].y += deltat * udashy; } time += deltat; } static void draw(){ Graphics2D g = (Graphics2D) img.getGraphics(); g.setColor(new Color(128, 212, 255)); g.fillRect(0, 0, 500, 400); int size = 8; double s = size / parcSize; for(int i = 0; i < elementCount; ++i){ Elem el = elements[i]; g.setColor(el.getColor()); g.fillOval((int)(el.x * s) + 30, (int)(el.y * s) + 20, size * 3/2, size * 3/2); } g.setColor(Color.BLACK); g.drawString(String.format("%2.4f", time), 10, 15); g.dispose(); label.repaint(); } }
