/** applet No. 1048 * * water h=1cm in SI unit - smoothed particle hydrodynamics 2D * - multi thread SPH2D :: display - asynchronous * * Created by Ikeuchi Mitsuru on October 28 2006. * Copyright (c) 2006-2007 Ikeuchi Mitsuru. All rights reserved. * * ver 0.0.1 2006.10.28 created * ver 0.0.2 2006.12.06 improved code * ver 0.0.3 2007.05.26 improved code * */ import java.awt.*; import java.awt.event.*; import java.applet.*; public class waterMtPSPH2D extends Applet implements MouseListener, MouseMotionListener, ItemListener, ActionListener, AdjustmentListener, Runnable { // ------------------------------------ preset field ----------- Thread th = null; // for run()-paint() thread SPH2D sph = null; DrawGraph2D dg2d = new DrawGraph2D(); // for event Button bt_reset, bt_startStop, bt_step, bt_go; Choice ch_view; Scrollbar sc_alpha, sc_dispScale; // for off-paint buffer Dimension dim; Image imgOff; Graphics gOff; int sleepTime = 100; int dispMode = 0; double viewScale = 1.0; int dgX, dgY, dgXb, dgYb; // mouse int thCount = 0; // ------------------------------ applet main thread ----------- public void init() { resize(630,420); setBackground(Color.white); dim = getSize(); imgOff = createImage(dim.width,dim.height); gOff = imgOff.getGraphics(); bt_reset= new Button("reset"); bt_reset.addActionListener(this); bt_startStop= new Button("start/stop"); bt_startStop.addActionListener(this); bt_step= new Button("step"); bt_step.addActionListener(this); bt_go = new Button("go"); bt_go.addActionListener(this); ch_view = new Choice(); ch_view.add("particle"); ch_view.add("velocity"); ch_view.add("density"); ch_view.add("temperature"); ch_view.add("pressure"); ch_view.addItemListener(this); ch_view.select("particle"); sc_alpha= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,510); sc_alpha.addAdjustmentListener(this); sc_dispScale= new Scrollbar(Scrollbar.HORIZONTAL,100,10,50,310); sc_dispScale.addAdjustmentListener(this); addMouseListener(this); addMouseMotionListener(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(bt_go); pnl.add(sc_alpha); pnl.add(sc_dispScale); pnl.add(ch_view); add(pnl,"North"); } public void start() { if (sph == null) { sph = new SPH2D(); sph.start(); } if (th == null) { th = new Thread(this); th.start(); } } public void stop() { if (th != null) th = null; if (sph != null) sph = null; } // ---------------------------------- event listener ----------- public void itemStateChanged(ItemEvent ev){ if (ev.getSource() == ch_view) { dispMode = ch_view.getSelectedIndex(); } } public void actionPerformed(ActionEvent ev){ if (ev.getSource() == bt_reset){ sph.reset(); thCount = 0; } else if (ev.getSource() == bt_startStop){ if (sph.getStartSW()==0) { sph.setStartSW(1); } else { sph.setStartSW(0); } } else if (ev.getSource() == bt_step){ if (sph.getStartSW()==0) sph.setStepSW(); } else if (ev.getSource() == bt_go){ sph.setFree(); } } public void adjustmentValueChanged(AdjustmentEvent ev){ if (ev.getSource() == sc_alpha) { sph.setAlpha( 0.00001*sc_alpha.getValue() ); } else if (ev.getSource() == sc_dispScale) { viewScale = 0.01*sc_dispScale.getValue(); dg2d.setViewScale(viewScale); } } public void mousePressed(MouseEvent ev){ } public void mouseReleased(MouseEvent ev){ dgXb = 0; dgYb = 0; dgX = 0; dgY = 0; } public void mouseClicked(MouseEvent ev){ } public void mouseEntered(MouseEvent ev){ } public void mouseExited (MouseEvent ev){ } public void mouseMoved (MouseEvent ev){ } public void mouseDragged(MouseEvent ev){ dgXb = dgX; dgYb = dgY; dgX=ev.getX(); dgY=ev.getY(); if (dgXb==0 && dgYb==0) { dgXb = dgX; dgYb = dgY; } dg2d.setShift( (dgX-dgXb), (dgY-dgYb) ); } // ========================= run() - paint() loop ============== public void run() { while (th != null) { try { Thread.sleep(sleepTime); } catch (InterruptedException e) { } thCount += 1; offPaint(); repaint(); } } public void update(Graphics g) { paint(g); } public void paint(Graphics g) { g.drawImage(imgOff,0,0,this); } void offPaint() { if (sph != null) { dg2d.plotGraph(gOff, sph, dispMode, thCount); } } } // =============================== draw graph 2D class =========== class DrawGraph2D { private Graphics gOff; private SPH2D sph; private int NNp; private double hh0 = SPH2D.hh0; private double xMax = SPH2D.xMax; private double yMax = SPH2D.yMax; private int gx0 = 30; private int gy0 = 80; private double xShift = 0.0; private double yShift = 0.0; private double dispScale = 400.0/xMax; private double viewScale = 1.0; private int xImageLoc = gx0+(int)(dispScale*viewScale*xShift); private int yImageLoc = gy0+(int)(dispScale*viewScale*yShift); private int drawnCount = 0; // ------------------------------- class constructor ----------- DrawGraph2D() { } // ----------------------------------- class methods ----------- public void plotGraph(Graphics gg, SPH2D sp, int dispMode, int thCount) { int sphCount; double sc; gOff = gg; sph = sp; if (sp==null) return; NNp = sph.getNNp(); sphCount = sph.getCount(); if (sphCount==drawnCount) return; drawnCount = sphCount; gOff.setColor(Color.white); gOff.fillRect(0,0,630,420); sc = dispScale*viewScale; gOff.setColor(Color.black); gOff.drawRect(xImageLoc,yImageLoc,(int)(sc*xMax+0.5),(int)(sc*yMax+0.5)); if (dispMode==0) { // oval plot ovalPlot(xImageLoc,yImageLoc,sc); } else if (dispMode==1) { // flow line linePlot(xImageLoc,yImageLoc,sc,50.0); } else if (dispMode==2) { // density fieldPlot(xImageLoc,yImageLoc,sc,0, 0.0002); } else if (dispMode==3) { // temperature fieldPlot(xImageLoc,yImageLoc,sc,1, 0.01); } else if (dispMode==4) { // pressure fieldPlot(xImageLoc,yImageLoc,sc,2, 0.0005); } gOff.setColor(Color.black); gOff.drawString("t="+(int)(sph.getTime()*1000*100.0+0.5)/100.0+" ms",630/6*0+10,40); if (sph.getStartSW()==1) { gOff.drawString("started",630/6*1+10,40); } else { gOff.drawString("stopped",630/6*1+10,40); } gOff.drawString("a="+(int)(sph.getAlpha()*100000.0+0.5)/100.0+"/1000",630/6*3+10,40); gOff.drawString("scale="+(int)(viewScale*100.0+0.5)+"%",630/6*4+10,40); gOff.drawString("view",630/6*5+10,40); gOff.drawString("maxSec="+sph.getMaxSection()+" ",10,60); gOff.drawString("Lx=60cm, Ly=40cm, h=1.0cm, water ",200,60); gOff.setColor(Color.black); gOff.drawString("thread MD = "+sph.getCount()+" ",20,400); gOff.drawString("thread disp = "+thCount+" ",220,400); } // ------------------------------------ plot methods ----------- void ovalPlot(int gx, int gy, double sc) { int i,j,ir,ix,iy,ip; ir = (int)(2.0*sph.gethh0()*sc+0.5); for (i=0; i=0) { if (sph.distance(i,ip)<0.5*xMax) { gOff.drawLine(gx+(int)(sph.getxx(i)*sc+0.5), gy+(int)(sph.getyy(i)*sc+0.5), gx+(int)(sph.getxx(ip)*sc+0.5), gy+(int)(sph.getyy(ip)*sc+0.5)); } } } } } } private void setColorOf(int i) { int k; k = sph.getKind(i); if (k>0) { gOff.setColor(Color.blue); } else if (k==0) { gOff.setColor(Color.green); } else { gOff.setColor(Color.orange); } } // void linePlot(int gx, int gy, double sc, double mag) { int i,ix,iy,ir; ir = (int)(2.0*sph.gethh0()*sc+0.5); for (i=0; i0.0) { gOff.setColor(Color.getHSBColor(0.66f,0.8f,0.8f)); } else { gOff.setColor(Color.getHSBColor(0.01f,0.8f,0.8f)); } gOff.drawLine(ix,iy,ix+(int)(sph.getvx(i)*viewScale*mag+0.5),iy+(int)(sph.getvy(i)*viewScale*mag+0.5)); if (sph.getKind(i)<=0) { setColorOf(i); gOff.drawOval(ix-ir/2,iy-ir/2, ir, ir); } } } // void fieldPlot(int gx, int gy, double sc, int mode, double mag) { int i,ir,ix,iy; double c; ir = (int)(2.0*sph.gethh0()*sc+0.5); for (i=0; i1.0) c = 1.0; gOff.setColor(Color.getHSBColor((float)(0.66-0.66*c),0.8f,0.8f)); gOff.drawOval(ix-ir/2,iy-ir/2, ir, ir); } } private double fieldColor(int i, int mode, double mag) { if (mode==0) { // density return mag*sph.getDensity(i); } else if (mode==1) { // temperature return mag*(sph.getTemperature(i)-273.0); } else if (mode==2) { // pressure return mag*sph.getPressure(i); } return 0.0; } // ------------------------------------- I/O methods ----------- public void setShift(double xDisplacement, double yDisplacement) { xShift += (1.0/viewScale/dispScale)*xDisplacement; yShift += (1.0/viewScale/dispScale)*yDisplacement; xImageLoc = gx0+(int)(dispScale*viewScale*xShift); yImageLoc = gy0+(int)(dispScale*viewScale*yShift); } public void setViewScale(double sc) { viewScale = sc; } public double getViewScale() { return viewScale; } } // =================== smoothed particle hydrodynamics =========== class SPH2D extends Thread { public static final double xMax = 0.6; // in m public static final double yMax = 0.4; // in m public static final double hh0 = 0.01; // in m private double t = 0.0; private double dt = 0.0005; // in s private double rBond = hh0*1.2; private double alpha = 0.001; private double thermalDiffusivity = 1.5e-7; // in m2/s private double dens0 = 1000.0; // in kg/m3 private int NNp = 5000; private int kind[] = new int [NNp]; private int bond[][] = new int[NNp][8]; private double xx[] = new double [NNp]; private double yy[] = new double [NNp]; private double vx[] = new double [NNp]; private double vy[] = new double [NNp]; private double ax[] = new double [NNp]; private double ay[] = new double [NNp]; private double hh[] = new double [NNp]; private double mass[] = new double [NNp]; private double cp[] = new double [NNp]; private double density[] = new double [NNp]; private double pressure[] = new double [NNp]; private double energy[] = new double [NNp]; private double power[] = new double [NNp]; private int Nsx = 60; private int Nsy = 40; private int section[][][] = new int[Nsx][Nsy][20]; private int resetSW = 0; private int startSW = 1; private int stepSW = 0; private int sleepTime = 1; private int stopSleepTime = 100; private int count = 0; // ------------------------------- class constructor ----------- SPH2D() { setInitialCondition(); } // -------------------------------------- thread run ----------- public void run() { while (true) { count += 1; timeEvolution(); if (sleepTime>0) { try { Thread.sleep(sleepTime); } catch (InterruptedException e) { } } } } // --------------------------- set initial condition ----------- private void setInitialCondition() { int i,j,ix,iy; double a; count = 0; t = 0.0; i = 0; // water; density=1.0g/cc, specific heat=4200J/kgK a = 0.011; for (ix=0; ix<30; ix++) { for (iy=0; iy<30; iy++) { kind[i] = 1; mass[i] = 0.15; cp[i] = 4200.0; xx[i] = 0.01+a*(ix)+0.5*a*(iy%2); yy[i] = 0.1+a*0.866*iy; vx[i] = 0.5*(Math.random()-0.5); vy[i] = 0.5*(Math.random()-0.5); ax[i] = 0.0; ay[i] = 0.0; hh[i] = hh0; i += 1; } } // solid body; density=0.5g/cc, specific heat=2000J/kgK for (ix=0; ix<10; ix++) { for (iy=0; iy<10; iy++) { kind[i] = 0; mass[i] = 0.075; cp[i] = 1000.0; xx[i] = 0.4+0.866*rBond*ix; yy[i] = 0.05+rBond*iy+0.5*rBond*(ix%2); vx[i] = 0.0; vy[i] = 0.0; ax[i] = 0.0; ay[i] = 0.0; i += 1; } } NNp = i; bonding(); dividingSection(); setDensity(); for (i=0; i0) { pressure[i] = 20.0*(Math.pow(density[i]/dens0,7.0)-1); // pressure[i] = 8.31*density[i]*(energy[i]/(mass[i]*cp[i])); } else { pressure[i] = 0.5*density[i]; } } for (i=0; i xMax) { xx[i] -= xMax; } if (yy[i] < 0.0) { yy[i] = 0.0; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; } if (yy[i] > yMax) { yy[i] = yMax; vx[i] = rr*vx[i]; vy[i] = -rr*vy[i]; } } } private void changeCalc() { int i,j,i0,i1,j0,j1,ii,jj,iq,iii,jjj; double ai,aj,b,xij,yij,r2,rij,gradW,gradWx,gradWy,lapW,a; double rcc,rcc2,vr,mu,cv,pij; cv = 100.0; rcc = hh0*2.0; rcc2 = rcc*rcc; for (i=0; i=Nsy) j1 = Nsy-1; for(ii=i0;ii<=i1;ii++) { iii = (ii+Nsx)%Nsx; for(jj=j0;jj<=j1;jj++) { jjj = jj; for(iq=1;iq<=section[iii][jjj][0];iq++) { j = section[iii][jjj][iq]; xij = xx[i] - xx[j]; yij = yy[i] - yy[j]; if (xij>0.5*xMax) xij -= xMax; if (xij<-0.5*xMax) xij += xMax; r2=xij*xij+yij*yij; if (r20) { if (kind[j]>0 && kind[i]>0) { aj = pressure[j]/(density[j]*density[j]); vr = xij*(vx[i]-vx[j])+yij*(vy[i]-vy[j]); pij = 0; if (vr<0.0) { mu = 0.5*vr/(r2+0.001*0.5*0.5); pij = (-alpha*cv*mu+0.0*mu*mu)/((density[i]+density[j])/2.0); } b = mass[j]*(ai+aj+pij); gradW = dwwvdr(rij, (hh[i]+hh[j])/2.0); gradWx = gradW*xij/rij; gradWy = gradW*yij/rij; ax[i] += -b*gradWx; ay[i] += -b*gradWy; lapW = d2wwvdr2(rij, (hh[i]+hh[j])/2.0); power[i] += 0.5*b*((vx[i]-vx[j])*gradWx+(vy[i]-vy[j])*gradWy) + thermalDiffusivity*mass[j]/density[j]*(energy[j]-energy[i])*lapW; } else { a = force(rij,i,j)/mass[i]; ax[i] += a*xij/rij; ay[i] += a*yij/rij; } } } } } } } } private double force(double r, int ii, int jj) { int i,j; double ret,rc,ra,ra1,ra2,ram,fa,f,r1,r3,kc; if (kind[ii]>=0 || kind[jj]>=0) { rc = hh0*1.2; ret = 0.0; if (rra2) { bondUnreg(i,j); bondUnreg(j,i); return ( 0.0 ); } f = 0.0; if (bondQ(i,j)==0) { if (r=ra2) { f = 0.0; bondUnreg(i,j); bondUnreg(j,i); } else if (ra1<=r && r xMax-rc) { r = xx[i]-(xMax-rc); ax[i] += -k*r*r; } */ if (yy[i] < rc) { r = rc-yy[i]; ay[i] += k*r*r; } if (yy[i] > yMax-rc) { r = yy[i]-(yMax-rc); ay[i] += -k*r*r; } } private void setDensity() { int i,j,i0,i1,j0,j1,ii,jj,iq,iii,jjj,ikind; double s,rij,r2,xij,yij; double rcc,rcc2; rcc = hh0*2.0; rcc2 = rcc*rcc; for (i=0; i=Nsy) j1 = Nsy-1; for(ii=i0;ii<=i1;ii++) { iii = (ii+Nsx)%Nsx; for(jj=j0;jj<=j1;jj++) { jjj = jj; for(iq=1;iq<=section[iii][jjj][0];iq++) { j = section[iii][jjj][iq]; if (kind[j]==ikind) { xij = xx[i] - xx[j]; yij = yy[i] - yy[j]; if (xij>0.5*xMax) xij -= xMax; if (xij<-0.5*xMax) xij += xMax; r2=xij*xij+yij*yij; if (r2=Nsx) i = Nsx-1; j = (int)(Nsy*yy[ip]/yMax); if (j>=Nsy) j = Nsy-1; iq = section[i][j][0]+1; section[i][j][0] = iq; section[i][j][iq] = ip; } } private int maxSection() { int i,j,m; m = 0; for(i=0;im) m = section[i][j][0]; } } return ( m ); } // smoothed particle private double kernel(double r, double h) { double a,q,ret; ret = 0.0; q = r/h; a = 10.0/(7.0*3.141592)/(h*h); if (q<1.0) { ret = a*(1.0-1.5*q*q+0.75*q*q*q); } else if (q<2.0) { ret = a*0.25*(2.0-q)*(2.0-q)*(2.0-q); } return ( ret ); } private double dwwvdr(double r, double h) { double a,q,ret; ret = 0.0; q = r/h; a = 10.0/(7.0*3.141592)/(h*h*h); if (q<1.0) { ret = a*(-3.0*q+2.25*q*q); } else if (q<=2.0) { ret = -a*0.75*(2.0-q)*(2.0-q); } return ( ret ); } private double d2wwvdr2(double r, double h) { double a,q,ret; ret = 0.0; q = r/h; a = 10.0/(7.0*3.141592)/(h*h*h*h); if (q<1.0) { ret = a*(-3.0+4.5*q); } else if (q<=2.0) { ret = a*0.75*(2.0-q); } return ( ret ); } // ----------------------------------------- utility ----------- // bonding manager private void bonding() { int i,j; double r2,rb,rb2; rb = rBond*1.2; rb2=rb*rb; bondInit(); for (i=0; irb2) { bondUnreg(ii,jj); bondUnreg(jj,ii); q = 0; } break; } } return (q); } private double rSq(int i,int j) { double xij,yij; xij = xx[i] - xx[j]; yij = yy[i] - yy[j]; if (xij>0.5*xMax) xij -= xMax; if (xij<-0.5*xMax) xij += xMax; return (xij*xij+yij*yij); } // object fixed(kind=0) -> move(kind=-1) private void gogogo() { int i; for (i=0; i