/**********************************************************************/ /* CFD.java */ /* A 2-D Inviscid Flow Simulation. */ /* Saleh Elmohamed. */ /**********************************************************************/ /* A note to experienced object oriented programmers: */ /* Some parts of this code will look sloppy because I used inline */ /* code for some things which could have been done more elegantly */ /* using public methods. */ /* This was done because inlining by hand increased the execution */ /* speed by about 20%. */ /**********************************************************************/ /*The simulation: *----------- * The code simulates a 2-D inviscid flow through an axisymmetric * nozzle. The simulation yields countour plots of all flow variables, * including velocity components, pressure, mach number, density and * entropy, and temperature. The plots show the location of any shock * wave that would reside in the nozzle. Also, the code finds the * steady state solution to the 2-D Euler equations * (see report and reference). * * It uses a 4-stage Runge-Kutta time-stepping algorithm and a finite * volume centeral-difference technique to find the solution. * At each stage, the residual is multiplied by a different value. * Generally, this gives a bit more accurate results since you go * over more than one stage. * A numerical dissipation model is emplyed to dampen any spurious * oscillations and prevent the soultion from blowing up in the * presence of shock waves (Jameson et al, 81). The dampening order * is set by: (1) compute the even derivatives (i.e. 2nd and 4th), * (2) set the order proportional to this. * * If the inlet and exit conditions are supersonic, the code requires * about 1800 iterations to complete execution. However, if subsonic * boundary conditions are employed, the number of iterations * required jumps to about 6000, as new errors are introduced into * the domain by these conditions (see Dang, 1997). * The flow is going from left to right across the nozzle. To see the * flow through, a "tuft" key is added as well as a "grid" key. * The flow can be either "subsonic", "transonic", or "supersonic". * The subsonic flow requires < 1, the transonic flow is partly sub * and partly super, and the supersonic requires a flow of greater * than and equal 1. * /* -------------------* / /* Known bug: */ /*******************************************************************/ /*** The temperature limits can show up as zero on the scale */ /* because the temperature is normalized to a small value. */ /* It's harmless, but irritating. */ /* * * References: ************ * 1. "Numerical Solutions of the Euler Equations by Finite Volume * Methods Using Runge-Kutta Time-Stepping Schemes", A. Jameson, * W. Schmidt, and E. Turkel, paper no: AIAA-81-1259, AIAA 14th * Fluid and Plasma Dynamics Conference, June 23-35, 1981, * Palo Alto, Calif. * * 2. "Monograph in CFD", T. Q. Dang, Dept. of Mechnical, Aerospace * and Manufacturing Engineering. Syracuse Univ. 1997. * ********************************************************************/ import java.applet.*; /* The usual stuff */ import java.awt.*; import java.io.InputStream; /* Lets you read files */ import java.io.StreamTokenizer; /* Lets you interprit what you read */ import java.net.*; /* Let you open a URL */ import java.io.IOException; /* Deals with io exceptions */ import java.io.*; import hpjava.hputil.*; import hpjava.adlib.*; import hpjava.lang.*; public class CFD extends HPApplet implements HPRunnable{ /* Declare "globals" */ Checkbox drawGridBox, pauseBox, tuftsBox; Choice plotChoice, timeChoice, viewChoice; Color ctable[]; DataPanel dataPanel; DampingPanel sdampingPanel, fdampingPanel; MachffPanel machffPanel; String inname = null; StringBuffer buffer = new StringBuffer(); Thread drawThread = null; double machffmin; double machffmax; final double dampingMin = 0.1; final double dampingMax = 5.0; Procs control; public void hpInit() { control = new Procs0(); InputStream instream = null; try { inname = "Tran.dat";//getParameter("inputfile"); instream = new URL(getDocumentBase(), inname).openStream(); dataPanel = new DataPanel(instream, control); } catch (IOException e) { System.err.println("Input/Output Exception: probably bad "+ "filename or bad connection." + e.getMessage()); } Panel controlPanelSouth = new Panel(); Panel controlPanelNorth = new Panel(); setLayout(new BorderLayout()); on(control){ /* Make the control panel */ controlPanelNorth.setLayout(new FlowLayout(FlowLayout.LEFT, 5, 5)); controlPanelSouth.setLayout(new FlowLayout(FlowLayout.LEFT, 5, 5)); plotChoice = new Choice(); plotChoice.addItem("Contours"); plotChoice.addItem("Line Plot"); viewChoice = new Choice(); viewChoice.addItem("Density"); viewChoice.addItem("Pressure"); viewChoice.addItem("Temperature"); viewChoice.addItem("Mach Number"); viewChoice.addItem("Axial Velocity"); viewChoice.addItem("Normal Velocity"); viewChoice.addItem("Total Pressure"); drawGridBox = new Checkbox("Grid", null, false); tuftsBox = new Checkbox("Tufts", null, false); timeChoice = new Choice(); timeChoice.addItem("Local Timestep"); timeChoice.addItem("Time Accurate"); pauseBox = new Checkbox("Pause", null, false); pauseBox.setFont(new Font("Helvetica", Font.PLAIN, 10)); drawGridBox.setFont(new Font("Helvetica", Font.PLAIN, 10)); tuftsBox.setFont(new Font("Helvetica", Font.PLAIN, 10)); } /* Set limits on mach slider. */ /* These shifting limits are necessary */ /* because there is a singularity */ /* in the boundary conditions when the */ /* mach number == 1.0 */ if (dataPanel.machff < 0.9) { machffmin = 0.2; machffmax = 0.9; } else if (dataPanel.machff > 1.3) { machffmin = 1.3; machffmax = 2.5; } else { machffmin = 0.2; machffmax = 0.9; } machffPanel = new MachffPanel("Inlet Mach Number", machffmin, machffmax, dataPanel.machff); sdampingPanel = new DampingPanel("2nd Order Damping", dampingMin, dampingMax, dataPanel.secondOrderDamping); fdampingPanel = new DampingPanel("4th Order Damping", dampingMin, dampingMax, dataPanel.fourthOrderDamping); dataPanel.setFont(new Font("Helvetica", Font.PLAIN, 10)); on(control) { /* Add the button panel */ controlPanelSouth.setFont(new Font("Helvetica", Font.PLAIN, 10)); Label label2 = new Label (" of"); controlPanelNorth.add(plotChoice); controlPanelNorth.add(label2); controlPanelNorth.add(viewChoice); controlPanelNorth.add(drawGridBox); controlPanelNorth.add(tuftsBox); controlPanelNorth.add(timeChoice); controlPanelSouth.add(machffPanel); controlPanelSouth.add(fdampingPanel); controlPanelSouth.add(sdampingPanel); add("South", controlPanelSouth); add("North", controlPanelNorth); } /* Add contour plot to main screen */ add("Center", dataPanel); } public void hpStart() { if (drawThread == null) { drawThread = HPutil.startThread(this, Thread.MIN_PRIORITY+2); } } public void hpRun(){ Adlib.init(); dataPanel.broadcastInitialData(); while (!Adlib.broadcast(shuttingDown,control)){ dataPanel.machff = Adlib.broadcast(machffPanel.sliderValue, control); dataPanel.secondOrderDamping = Adlib.broadcast(sdampingPanel.sliderValue, control); dataPanel.fourthOrderDamping = Adlib.broadcast(fdampingPanel.sliderValue, control); dataPanel.nplot = Adlib.broadcast(nplot, control); dataPanel.nview = Adlib.broadcast(nview, control); dataPanel.ntime = Adlib.broadcast(ntime, control); dataPanel.setGridSwitch(Adlib.broadcast(gridState, control)); dataPanel.setTuftsSwitch(Adlib.broadcast(tuftState, control)); dataPanel.doIteration(); } Adlib.finish(); synchronized(this) { notify(); } } boolean shuttingDown = false; public synchronized void hpStop() { on(control){ shuttingDown = true; } try { wait(); }catch (InterruptedException e) { System.err.println("Execution interrupted by another thread."); System.err.println(e.getMessage()); } } int nplot, nview, ntime; boolean gridState, tuftState, pause ; public boolean hpHandleEvent(Event e) { if (e.target instanceof Choice) { nplot = plotChoice.getSelectedIndex(); nview = viewChoice.getSelectedIndex(); ntime = timeChoice.getSelectedIndex(); return true; } else if (e.target instanceof Checkbox) { gridState = drawGridBox.getState(); tuftState = tuftsBox.getState(); return true; } return false; } } /* The fundamental class is the data */ /* This executes when you create a new object */ class datafield implements HPspmd{ private double[[-,-]] a; /* Grid cell area */ double[[-,-]] deltat; /* Timestep */ private double machff; /* Inflow mach number */ public double machfftarget; /* Target inflow mach number */ public double secondOrderDamping; public double fourthOrderDamping; public int ntime; private double[[-,-]] opg, pg, pg1; /* Pressure */ private double[[-,-]] sxi, seta; private double[[-,-]] tg, tg1; /* Temperature */ private double[[-,-]] xnode, ynode; /* Storage of node coordinates */ double xmin=Double.MAX_VALUE, xmax=Double.MIN_VALUE, ymin=Double.MAX_VALUE, ymax=Double.MIN_VALUE; /* Max and Min of node coordinates */ double cff, uff, vff, pff, rhoff, tff, jplusff, jminusff; /* Far field values */ double datamax, datamin; int iter; int imax, jmax; /* Number of nodes in x and y direction */ private Statevector[[-,-]] d; /* Damping coefficients */ private Statevector[[-,-]] f, g; /* Flux Vectors */ private Statevector[[-,-]] r, ug1; private Statevector[[-,-]] ug; /* Storage of data */ private Statevector[[-,-]] tempug; final double Cp = 1004.5; /* specific heat, const pres. */ final double Cv=717.5; /* specific heat, const vol. */ final double gamma=1.4; /* Ratio of specific heats */ final double rgas=287.0; /* Gas Constant */ final int keytextwidth=50; /* Width of text field in color table */ final int ncontour = 16; /* Number of contours between min and max */ final double fourthOrderNormalizer = 0.02; /* Damping coefficients */ final double secondOrderNormalizer = 0.02; private double[[*]] xnodeE, xnodeW, xnodeN, xnodeS, ynodeE, ynodeW, ynodeN,ynodeS; Range x; Range y; Procs control; Procs p; double[[*,*]] seqXnode; double[[*,*]] seqYnode; Statevector[[*,*]] seqUg; double[[*,*]] seqTg; double[[*,*]] seqPg; datafield(InputStream instream, Procs control) throws IOException { this.control = control; p = new Procs2(2,2); int i, j, n; /* Dummy counters */ double scrap = 0.0; double scrap2 = 0; /* Temporary storage */ /* Convert the stream into tokens (which helps you parse it) *-----------------*/ // StreamTokenizer intokens = new StreamTokenizer(instream); Reader reader; StreamTokenizer intokens = null; // read header from the root thread. reader = new BufferedReader(new InputStreamReader(instream)); intokens = new StreamTokenizer(reader); /* Read header */ if (intokens.nextToken() == StreamTokenizer.TT_NUMBER) imax = (int) intokens.nval; else throw new IOException(); if (intokens.nextToken() == StreamTokenizer.TT_NUMBER) jmax = (int) intokens.nval; else throw new IOException(); if (intokens.nextToken() == StreamTokenizer.TT_NUMBER) { machfftarget = intokens.nval; machff = machfftarget; } else throw new IOException(); if (intokens.nextToken() == StreamTokenizer.TT_NUMBER) iter = (int) intokens.nval; else throw new IOException(); // read values to the sequantial array on the root. seqXnode = new double[[imax+2, jmax+2]] on control; seqYnode = new double[[imax+2, jmax+2]] on control; seqUg = new Statevector[[imax+2, jmax+2]] on control; seqTg = new double[[imax+2, jmax+2]] on control; seqPg = new double[[imax+2, jmax+2]] on control; on(control) { for(i = 0; i < imax+2; i++) for(j = 0; j < jmax+2; j++) seqUg[i,j] = new Statevector(); /* Now load the data */ for (i = 0; i < imax; ++i) for (j = 0; j < jmax; ++j) { if (intokens.nextToken() == StreamTokenizer.TT_NUMBER) { seqXnode[i,j] = (double) intokens.nval; xmax = Math.max(xmax, seqXnode[i,j]); xmin = Math.min(xmin, seqXnode[i,j]); if (intokens.nextToken() == StreamTokenizer.TT_NUMBER) { seqYnode[i,j] = (double) intokens.nval; ymax = Math.max(ymax, seqYnode[i,j]); ymin = Math.min(ymin, seqYnode[i,j]); if(intokens.nextToken() == StreamTokenizer.TT_NUMBER) { seqUg[i,j].a =(double) intokens.nval; if(intokens.nextToken() == StreamTokenizer.TT_NUMBER) { seqUg[i,j].b =(double) intokens.nval; if(intokens.nextToken() == StreamTokenizer.TT_NUMBER) { seqUg[i,j].c =(double) intokens.nval; if(intokens.nextToken() == StreamTokenizer.TT_NUMBER) { seqUg[i,j].d =(double)intokens.nval; scrap = seqUg[i,j].c/seqUg[i,j].a; scrap2 = seqUg[i,j].b/seqUg[i,j].a; seqTg[i,j] = seqUg[i,j].d/seqUg[i,j].a - (0.5 * (scrap * scrap + scrap2 * scrap2)); seqTg[i,j] = seqTg[i,j] / Cv; seqPg[i,j] = rgas * seqUg[i,j].a * seqTg[i,j]; } } } } } }else throw new IOException(); } } } public void broadcastInitialData(){ x = new ExtBlockRange(imax+2, p.dim(0), 2, 2); y = new ExtBlockRange(jmax+2, p.dim(1), 2, 2); deltat = new double [[x,y]] on p; opg = new double [[x,y]] on p; pg = new double [[x,y]] on p; pg1 = new double [[x,y]] on p; sxi = new double [[x,y]] on p; seta = new double [[x,y]] on p; tg = new double [[x,y]] on p; tg1 = new double [[x,y]] on p; ug = new Statevector[[x,y]] on p; tempug = new Statevector[[x,y]] on p; a = new double [[x,y]] on p; d = new Statevector[[x,y]] on p; f = new Statevector[[x,y]] on p; g = new Statevector[[x,y]] on p; r = new Statevector[[x,y]] on p; ug1 = new Statevector[[x,y]] on p; xnode = new double [[x,y]] on p; ynode = new double [[x,y]] on p; xnodeE = new double[[imax+2]] on p; xnodeW = new double[[imax+2]] on p; xnodeN = new double[[jmax+2]] on p; xnodeS = new double[[jmax+2]] on p; ynodeE = new double[[imax+2]] on p; ynodeW = new double[[imax+2]] on p; ynodeN = new double[[jmax+2]] on p; ynodeS = new double[[jmax+2]] on p; on(p) { overall(i = x for :) overall(j = y for :){ d[i,j] = new Statevector(); f[i,j] = new Statevector(); g[i,j] = new Statevector(); r[i,j] = new Statevector(); ug[i,j] = new Statevector(); ug1[i,j] = new Statevector(); } } /* Set farfield values (we use normalized units for everything */ cff = 1.0; vff = 0.0; pff = 1.0 / gamma; rhoff = 1.0; tff = pff / (rhoff * rgas); // remap sequantial array which read by root theard to distributed array. Adlib.remap(xnode, seqXnode); Adlib.remap(ynode, seqYnode); Adlib.remap(xnodeE, xnode[[:, 0]]); Adlib.remap(xnodeW, xnode[[:, jmax-1]]); Adlib.remap(xnodeN, xnode[[0, :]]); Adlib.remap(xnodeS, xnode[[imax-1, :]]); Adlib.remap(ynodeE, ynode[[:, 0]]); Adlib.remap(ynodeW, ynode[[:, jmax-1]]); Adlib.remap(ynodeN, ynode[[0, :]]); Adlib.remap(ynodeS, ynode[[imax-1, :]]); Adlib.remap(ug[[1:imax, 1:jmax]], seqUg[[0:imax-1, 0:jmax-1]]); Adlib.remap(tg[[1:imax, 1:jmax]], seqTg[[0:imax-1, 0:jmax-1]]); Adlib.remap(pg[[1:imax, 1:jmax]], seqPg[[0:imax-1, 0:jmax-1]]); // breadcast values read by root thread. xmax = Adlib.broadcast(xmax, control); xmin = Adlib.broadcast(xmin, control); ymax = Adlib.broadcast(ymax, control); ymin = Adlib.broadcast(ymin, control); /* Calculate grid cell areas */ Adlib.writeHalo(xnode, wlo, whi); Adlib.writeHalo(ynode, wlo, whi); on(p) { overall (i = x for 1 : imax-1) overall (j = y for 1 : jmax-1) { a[i,j] = 0.5 * ((xnode[i,j] - xnode[i-1,j-1]) * (ynode[i-1,j] - ynode[i,j-1])- (ynode[i,j] - ynode[i-1,j-1]) * (xnode[i-1,j] - xnode[i,j-1])); } } } int wlo[] = {2, 2}, whi[] = {2, 2}; public void drawKey(Color ctable[], Graphics g, int xkeywidth, int ykeywidth, int xcanvaswidth, int ycanvaswidth) { on(control) { int colorbandheight, colorbandwidth; int n; int yband; int xoffset, yoffset; FontMetrics keyFontMetrics; xoffset = xcanvaswidth - xkeywidth; yoffset = ycanvaswidth - ykeywidth; colorbandwidth = xcanvaswidth - keytextwidth - (xoffset+5); colorbandheight = (ykeywidth-10) / (ctable.length-1) + 1; g.setColor(Color.black); g.fillRect(xoffset, yoffset, xkeywidth, ykeywidth); /* Draw the colors themselves */ for (n = 0; n < ctable.length; ++n) { g.setColor(ctable[n]); yband = ycanvaswidth - 5 - (int) (((double)(n+1)/(double) ctable.length) * (double)(ykeywidth-10)); g.fillRect(xoffset+5, yband, colorbandwidth, colorbandheight); } /* Write in the text for the key */ keyFontMetrics = g.getFontMetrics(); keyFontMetrics.getHeight(); g.setColor(Color.white); StringBuffer writeBuffer = new StringBuffer().append(datamin); writeBuffer.setLength(4); g.drawString(" "+writeBuffer, xoffset+5+colorbandwidth, ycanvaswidth - 5); StringBuffer writeBuffer2 = new StringBuffer().append(datamax); writeBuffer2.setLength(4); g.drawString(" "+writeBuffer2, xoffset+5+colorbandwidth, keyFontMetrics.getHeight() + 5); /* Draw outline frame */ g.setColor(Color.white); g.drawRect(xoffset, yoffset, xkeywidth, ykeywidth); } } void doIteration() { on(p) { double error; double scrap; int i, j; /* Record the old pressure values */ HPutil.copy(opg[[1:imax-1, 1:jmax-1]], pg[[1:imax-1, 1:jmax-1]]); /* Set new far field mach number (if necessary) */ if (machff != machfftarget) { if (Math.abs(machfftarget-machff) <= 0.02) machff = machfftarget; else { if (machfftarget > machff) { machff += 0.02; } else machff -= 0.02; } } calculateDummyCells(pg, tg, ug); calculateDeltaT(); calculateDamping(pg, ug); /* Do the integration */ /*...... Step 1 ......*/ calculateF(pg, tg, ug); calculateG(pg, tg, ug); calculateR(); overall(i = x for 1 : imax - 1) overall(j = y for 1 : jmax - 1) { ug1[i,j].a=ug[i,j].a-0.25*deltat[i,j]/a[i,j]* (r[i,j].a-d[i,j].a); ug1[i,j].b=ug[i,j].b-0.25*deltat[i,j]/a[i,j]* (r[i,j].b-d[i,j].b); ug1[i,j].c=ug[i,j].c-0.25*deltat[i,j]/a[i,j]* (r[i,j].c-d[i,j].c); ug1[i,j].d=ug[i,j].d-0.25*deltat[i,j]/a[i,j]* (r[i,j].d-d[i,j].d); } calculateStateVar(pg1, tg1, ug1); /*...... Step 2 ......*/ calculateDummyCells(pg1, tg1, ug1); calculateF(pg1, tg1, ug1); calculateG(pg1, tg1, ug1); calculateR(); overall(i = x for 1 : imax - 1) overall(j = y for 1 : jmax - 1) { ug1[i,j].a= ug[i,j].a-0.33333*deltat[i,j]/a[i,j]*(r[i,j].a-d[i,j].a); ug1[i,j].b= ug[i,j].b-0.33333*deltat[i,j]/a[i,j]*(r[i,j].b-d[i,j].b); ug1[i,j].c= ug[i,j].c-0.33333*deltat[i,j]/a[i,j]*(r[i,j].c-d[i,j].c); ug1[i,j].d= ug[i,j].d-0.33333*deltat[i,j]/a[i,j]*(r[i,j].d-d[i,j].d); } calculateStateVar(pg1, tg1, ug1); /*...... Step 3 ......*/ calculateDummyCells(pg1, tg1, ug1); calculateF(pg1, tg1, ug1); calculateG(pg1, tg1, ug1); calculateR(); overall(i = x for 1 : imax - 1) overall(j = y for 1 : jmax - 1) { ug1[i,j].a= ug[i,j].a-0.5*deltat[i,j]/a[i,j]*(r[i,j].a-d[i,j].a); ug1[i,j].b= ug[i,j].b-0.5*deltat[i,j]/a[i,j]*(r[i,j].b-d[i,j].b); ug1[i,j].c= ug[i,j].c-0.5*deltat[i,j]/a[i,j]*(r[i,j].c-d[i,j].c); ug1[i,j].d= ug[i,j].d-0.5*deltat[i,j]/a[i,j]*(r[i,j].d-d[i,j].d); } calculateStateVar(pg1, tg1, ug1); /*..... Step 4 (final step) .....*/ calculateDummyCells(pg1, tg1, ug1); calculateF(pg1, tg1, ug1); calculateG(pg1, tg1, ug1); calculateR(); overall(i = x for 1 : imax - 1) overall(j = y for 1 : jmax - 1) { ug[i,j].a -= deltat[i,j]/a[i,j]*(r[i,j].a-d[i,j].a); ug[i,j].b -= deltat[i,j]/a[i,j]*(r[i,j].b-d[i,j].b); ug[i,j].c -= deltat[i,j]/a[i,j]*(r[i,j].c-d[i,j].c); ug[i,j].d -= deltat[i,j]/a[i,j]*(r[i,j].d-d[i,j].d); } calculateStateVar(pg, tg, ug); /* calculate RMS Pressure Error. */ error = 0.0; overall(i = x for 1 : imax - 1) overall(j = y for 1 : jmax - 1) { scrap = pg[i,j]-opg[i,j]; error += scrap*scrap; } error = Math.sqrt(error / (double)((imax-1) * (jmax-1)) ); } } private void calculateStateVar(double [[-,-]] localpg, double [[-,-]] localtg, Statevector[[-,-]] localug) { /* Calculates the new state values for Range-Kutta */ /* Works for default values, 8/15 at 9:45 pm */ double temp, temp2; int i, j; overall(i = x for 1 : imax - 1){ overall(j = y for 1 : jmax - 1) { temp = localug[i,j].b; temp2 = localug[i,j].c; localtg[i,j] = localug[i,j].d/localug[i,j].a - 0.5 * (temp*temp + temp2*temp2)/(localug[i,j].a*localug[i,j].a); localtg[i,j] = localtg[i,j] / Cv; localpg[i,j] = localug[i,j].a * rgas * localtg[i,j]; } } } private void calculateR() { /* Works for default values, straight channel (all 0's) 8/11, 9:15 pm */ double deltax, deltay; double temp; int i,j; Statevector scrap; Adlib.writeHalo(f, wlo, whi); Adlib.writeHalo(g, wlo, whi); overall (i = x for 1 : imax-1){ overall (j = y for 1 : jmax-1) { /* Start by clearing R */ r[i,j].a = 0.0; r[i,j].b = 0.0; r[i,j].c = 0.0; r[i,j].d = 0.0; /* East Face */ deltay = (ynode[i,j] - ynode[i,j-1]); deltax = (xnode[i,j] - xnode[i,j-1]); temp = 0.5 * deltay; r[i,j].a += temp*(f[i,j].a + f[i+1,j].a); r[i,j].b += temp*(f[i,j].b + f[i+1,j].b); r[i,j].c += temp*(f[i,j].c + f[i+1,j].c); r[i,j].d += temp*(f[i,j].d + f[i+1,j].d); temp = -0.5*deltax; r[i,j].a += temp * (g[i,j].a+g[i+1,j].a); r[i,j].b += temp * (g[i,j].b+g[i+1,j].b); r[i,j].c += temp * (g[i,j].c+g[i+1,j].c); r[i,j].d += temp * (g[i,j].d+g[i+1,j].d); /* South Face */ deltay = (ynode[i,j-1] - ynode[i-1,j-1]); deltax = (xnode[i,j-1] - xnode[i-1,j-1]); temp = 0.5 * deltay; r[i,j].a += temp*(f[i,j].a+f[i,j-1].a); r[i,j].b += temp*(f[i,j].b+f[i,j-1].b); r[i,j].c += temp*(f[i,j].c+f[i,j-1].c); r[i,j].d += temp*(f[i,j].d+f[i,j-1].d); temp = -0.5*deltax; r[i,j].a += temp * (g[i,j].a+g[i,j-1].a); r[i,j].b += temp * (g[i,j].b+g[i,j-1].b); r[i,j].c += temp * (g[i,j].c+g[i,j-1].c); r[i,j].d += temp * (g[i,j].d+g[i,j-1].d); /* West Face */ deltay = (ynode[i-1,j-1] - ynode[i-1,j]); deltax = (xnode[i-1,j-1] - xnode[i-1,j]); temp = 0.5 * deltay; r[i,j].a += temp*(f[i,j].a+f[i-1,j].a); r[i,j].b += temp*(f[i,j].b+f[i-1,j].b); r[i,j].c += temp*(f[i,j].c+f[i-1,j].c); r[i,j].d += temp*(f[i,j].d+f[i-1,j].d); temp = -0.5*deltax; r[i,j].a += temp * (g[i,j].a+g[i-1,j].a); r[i,j].b += temp * (g[i,j].b+g[i-1,j].b); r[i,j].c += temp * (g[i,j].c+g[i-1,j].c); r[i,j].d += temp * (g[i,j].d+g[i-1,j].d); /* North Face */ deltay = (ynode[i-1,j] - ynode[i,j]); deltax = (xnode[i-1,j] - xnode[i,j]); temp = 0.5 * deltay; r[i,j].a += temp*(f[i,j].a+f[i+1,j].a); r[i,j].b += temp*(f[i,j].b+f[i+1,j].b); r[i,j].c += temp*(f[i,j].c+f[i+1,j].c); r[i,j].d += temp*(f[i,j].d+f[i+1,j].d); temp = -0.5*deltax; r[i,j].a += temp * (g[i,j].a+g[i,j+1].a); r[i,j].b += temp * (g[i,j].b+g[i,j+1].b); r[i,j].c += temp * (g[i,j].c+g[i,j+1].c); r[i,j].d += temp * (g[i,j].d+g[i,j+1].d); } } } private void calculateG(double[[-,-]] localpg, double[[-,-]] localtg, Statevector[[-,-]] localug) { /* Works for default values 8/15, 5:15 pm */ double temp, temp2, temp3; double v; int i, j; overall(i = x for 0 : imax) { overall(j = y for 0 : jmax) { v = localug[i,j].c / localug[i,j].a; g[i,j].a = localug[i,j].c; g[i,j].b = localug[i,j].b * v; g[i,j].c = localug[i,j].c*v + localpg[i,j]; temp = localug[i,j].b * localug[i,j].b; temp2 = localug[i,j].c * localug[i,j].c; temp3 = localug[i,j].a * localug[i,j].a; g[i,j].d = localug[i,j].c * (Cp * localtg[i,j]+ (0.5 * (temp + temp2)/(temp3))); } } } private void calculateF(double[[-,-]] localpg, double[[-,-]] localtg, Statevector[[-,-]] localug) { /* Works for default values 8/15, 4:50 pm */ double u; double temp1, temp2, temp3; int i, j; overall(i = x for 0 : imax) { overall(j = y for 0 : jmax) { u = localug[i,j].b/ localug[i,j].a; f[i,j].a = localug[i,j].b; f[i,j].b = localug[i,j].b *u + localpg[i,j]; f[i,j].c = localug[i,j].c * u; temp1 = localug[i,j].b * localug[i,j].b; temp2 = localug[i,j].c * localug[i,j].c; temp3 = localug[i,j].a * localug[i,j].a; f[i,j].d = localug[i,j].b * (Cp * localtg[i,j] + (0.5 * (temp1 + temp2)/(temp3))); } } } private void calculateDamping(double[[-,-]] localpg, Statevector[[-,-]] localug) { double adt, sbar; double nu2; double nu4; double tempdouble; int ascrap, i, j; Statevector temp = new Statevector(); Statevector temp2 = new Statevector(); Statevector scrap2 = new Statevector(), scrap4 = new Statevector(); nu2 = secondOrderDamping * secondOrderNormalizer; nu4 = fourthOrderDamping * fourthOrderNormalizer; /* First do the pressure switches */ /* Checked and works with defaults. */ Adlib.writeHalo(localpg, wlo, whi); overall (i = x for 1 : imax-1) overall (j = y for 1 : jmax-1) { sxi[i,j] = Math.abs(localpg[i+1,j] - 2.0 * localpg[i,j] + localpg[i-1,j]) / localpg[i,j]; seta[i,j] = Math.abs(localpg[i,j+1] - 2.0 * localpg[i,j] + localpg[i,j-1]) / localpg[i,j]; } /* Then calculate the fluxes ..... */ Adlib.writeHalo(localug, wlo, whi); Adlib.writeHalo(a, wlo, whi); Adlib.writeHalo(deltat, wlo, whi); Adlib.writeHalo(sxi, wlo, whi); Adlib.writeHalo(seta, wlo, whi); overall (i = x for 1 : imax-1){ overall (j = y for 1 : jmax-1) { /* Clear values */ /* East Face */ if (i` > 1 && i` < imax-1) { adt = (a[i,j] + a[i+1,j]) / (deltat[i,j] + deltat[i+1,j]); sbar = (sxi[i+1,j] + sxi[i,j]) * 0.5; } else { adt = a[i,j]/deltat[i,j]; sbar = sxi[i,j]; } tempdouble = nu2*sbar*adt; scrap2.a = tempdouble * (localug[i+1,j].a-localug[i,j].a); scrap2.b = tempdouble * (localug[i+1,j].b-localug[i,j].b); scrap2.c = tempdouble * (localug[i+1,j].c-localug[i,j].c); scrap2.d = tempdouble * (localug[i+1,j].d-localug[i,j].d); if (i` > 1 && i` < imax-1) { temp = localug[i+2,j].svect(localug[i-1,j]); temp2.a = 3.0*(localug[i,j].a-localug[i+1,j].a); temp2.b = 3.0*(localug[i,j].b-localug[i+1,j].b); temp2.c = 3.0*(localug[i,j].c-localug[i+1,j].c); temp2.d = 3.0*(localug[i,j].d-localug[i+1,j].d); tempdouble = -nu4*adt; scrap4.a = tempdouble*(temp.a+temp2.a); scrap4.b = tempdouble*(temp.a+temp2.b); scrap4.c = tempdouble*(temp.a+temp2.c); scrap4.d = tempdouble*(temp.a+temp2.d); } else { scrap4.a = 0.0; scrap4.b = 0.0; scrap4.c = 0.0; scrap4.d = 0.0; } temp.a = scrap2.a + scrap4.a; temp.b = scrap2.b + scrap4.b; temp.c = scrap2.c + scrap4.c; temp.d = scrap2.d + scrap4.d; d[i,j] = temp; /* West Face */ if(i` > 1 && i` < imax-1) { adt = (a[i,j] + a[i-1,j]) / (deltat[i,j] + deltat[i-1,j]); sbar = (sxi[i,j] + sxi[i-1,j]) *0.5; } else { adt = a[i,j]/deltat[i,j]; sbar = sxi[i,j]; } tempdouble = -nu2*sbar*adt; scrap2.a = tempdouble * (localug[i,j].a-localug[i-1,j].a); scrap2.b = tempdouble * (localug[i,j].b-localug[i-1,j].b); scrap2.c = tempdouble * (localug[i,j].c-localug[i-1,j].c); scrap2.d = tempdouble * (localug[i,j].d-localug[i-1,j].d); if (i` > 1 && i` < imax-1) { temp = localug[i+1,j].svect(localug[i-2,j]); temp2.a = 3.0*(localug[i-1,j].a-localug[i,j].a); temp2.b = 3.0*(localug[i-1,j].b-localug[i,j].b); temp2.c = 3.0*(localug[i-1,j].c-localug[i,j].c); temp2.d = 3.0*(localug[i-1,j].d-localug[i,j].d); tempdouble = nu4*adt; scrap4.a = tempdouble*(temp.a+temp2.a); scrap4.b = tempdouble*(temp.a+temp2.b); scrap4.c = tempdouble*(temp.a+temp2.c); scrap4.d = tempdouble*(temp.a+temp2.d); } else { scrap4.a = 0.0; scrap4.b = 0.0; scrap4.c = 0.0; scrap4.d = 0.0; } d[i,j].a += scrap2.a + scrap4.a; d[i,j].b += scrap2.b + scrap4.b; d[i,j].c += scrap2.c + scrap4.c; d[i,j].d += scrap2.d + scrap4.d; /* North Face */ if (j` > 1 && j` < jmax-1) { adt = (a[i,j] + a[i,j+1]) / (deltat[i,j] + deltat[i,j+1]); sbar = (seta[i,j] + seta[i,j+1]) * 0.5; } else { adt = a[i,j]/deltat[i,j]; sbar = seta[i,j]; } tempdouble = nu2*sbar*adt; scrap2.a = tempdouble * (localug[i,j+1].a-localug[i,j].a); scrap2.b = tempdouble * (localug[i,j+1].b-localug[i,j].b); scrap2.c = tempdouble * (localug[i,j+1].c-localug[i,j].c); scrap2.d = tempdouble * (localug[i,j+1].d-localug[i,j].d); if (j` > 1 && j` < jmax-1) { temp = localug[i,j+2].svect(localug[i,j-1]); temp2.a = 3.0*(localug[i,j].a-localug[i,j+1].a); temp2.b = 3.0*(localug[i,j].b-localug[i,j+1].b); temp2.c = 3.0*(localug[i,j].c-localug[i,j+1].c); temp2.d = 3.0*(localug[i,j].d-localug[i,j+1].d); tempdouble = -nu4*adt; scrap4.a = tempdouble*(temp.a+temp2.a); scrap4.b = tempdouble*(temp.a+temp2.b); scrap4.c = tempdouble*(temp.a+temp2.c); scrap4.d = tempdouble*(temp.a+temp2.d); } else { scrap4.a = 0.0; scrap4.b = 0.0; scrap4.c = 0.0; scrap4.d = 0.0; } d[i,j].a += scrap2.a + scrap4.a; d[i,j].b += scrap2.b + scrap4.b; d[i,j].c += scrap2.c + scrap4.c; d[i,j].d += scrap2.d + scrap4.d; /* South Face */ if (j` > 1 && j` < jmax-1) { adt = (a[i,j] + a[i,j-1]) / (deltat[i,j] + deltat[i,j-1]); sbar = (seta[i,j] + seta[i,j-1]) * 0.5; } else { adt = a[i,j]/deltat[i,j]; sbar = seta[i,j]; } tempdouble = -nu2*sbar*adt; scrap2.a = tempdouble * (localug[i,j].a-localug[i,j-1].a); scrap2.b = tempdouble * (localug[i,j].b-localug[i,j-1].b); scrap2.c = tempdouble * (localug[i,j].c-localug[i,j-1].c); scrap2.d = tempdouble * (localug[i,j].d-localug[i,j-1].d); if (j` > 1 && j` < jmax-1) { temp = localug[i,j+1].svect(localug[i,j-2]); temp2.a = 3.0*(localug[i,j-1].a-localug[i,j].a); temp2.b = 3.0*(localug[i,j-1].b-localug[i,j].b); temp2.c = 3.0*(localug[i,j-1].c-localug[i,j].c); temp2.d = 3.0*(localug[i,j-1].d-localug[i,j].d); tempdouble = nu4*adt; scrap4.a = tempdouble*(temp.a+temp2.a); scrap4.b = tempdouble*(temp.a+temp2.b); scrap4.c = tempdouble*(temp.a+temp2.c); scrap4.d = tempdouble*(temp.a+temp2.d); } else { scrap4.a = 0.0; scrap4.b = 0.0; scrap4.c = 0.0; scrap4.d = 0.0; } d[i,j].a += scrap2.a + scrap4.a; d[i,j].b += scrap2.b + scrap4.b; d[i,j].c += scrap2.c + scrap4.c; d[i,j].d += scrap2.d + scrap4.d; } } } private void calculateDeltaT() { double xeta, yeta, xxi, yxi; /* Local change in x and y */ int i, j; double mint; double c, q, r; double safety_factor = 0.3; overall(i = x for 1 : imax - 1) overall(j = y for 1 : jmax - 1){ xxi = (xnode[i,j] - xnode[i-1,j] + xnode[i,j-1] - xnode[i-1,j-1]) * 0.5; yxi = (ynode[i,j] - ynode[i-1,j] + ynode[i,j-1] - ynode[i-1,j-1]) * 0.5; xeta = (xnode[i,j] - xnode[i,j-1] + xnode[i-1,j] - xnode[i-1,j-1]) * 0.5; yeta = (ynode[i,j] - ynode[i,j-1] + ynode[i-1,j] - ynode[i-1,j-1]) * 0.5; q = (yeta * ug[i,j].b - xeta * ug[i,j].c); r = (-yxi * ug[i,j].b + xxi * ug[i,j].c); c = Math.sqrt (gamma * rgas * tg[i,j]); deltat[i,j] = safety_factor * 2.8284 * a[i,j] / ( (Math.abs(q) + Math.abs(r))/ug[i,j].a + c * Math.sqrt(xxi*xxi + yxi*yxi + xeta*xeta + yeta*yeta + 2.0 * Math.abs(xeta*xxi + yeta*yxi))); } /* If that's the user's choice, make it time accurate */ if (ntime == 1) HPutil.init(deltat[[1:imax-1, 1:jmax-1]], Adlib.minval(deltat)); } private void calculateDummyCells(double[[-,-]] localpg, double[[-,-]] localtg, Statevector[[-,-]] localug) { double c = 0; double jminus=0; double jplus=0; double s=0; double rho=0; double temp = 0; double u = 0; double v = 0; double scrap = 0; double scrap2= 0; double theta = 0; double uprime = 0; Vector2 norm = new Vector2(); Vector2 tan = new Vector2(); Vector2 u1 = new Vector2(); uff = machff; jplusff = uff + 2.0 / (gamma - 1.0) * cff; jminusff = uff - 2.0 / (gamma - 1.0) * cff; Adlib.writeHalo(localtg, wlo, whi); Adlib.writeHalo(localug, wlo, whi); Adlib.writeHalo(localpg, wlo, whi); overall(i = x for 1 : imax - 1) { at(j = y[0]){ /* Bottom wall boundary cells */ /* Routine checked by brute force for initial conditions. */ /* Routine checked by brute force for random conditions. */ /* Construct tangent vectors */ tan.ihat = xnode[i,j] - xnode[i-1,j]; tan.jhat = ynode[i,j] - ynode[i-1,j]; norm.ihat = - (ynode[i,j] - ynode[i-1,j]); norm.jhat = xnode[i,j] - xnode[i-1,j]; } scrap = tan.magnitude(); tan.ihat = tan.ihat / scrap; tan.jhat = tan.jhat / scrap; scrap = norm.magnitude(); norm.ihat = norm.ihat / scrap; norm.jhat = norm.jhat / scrap; at(j = y[0]) { /* Now set some state variables */ rho = localug[i,j+1].a; localtg[i,j] = localtg[i,j+1]; u1.ihat = localug[i,j+1].b / rho; u1.jhat = localug[i,j+1].c / rho; } u = u1.dot(tan) + u1.dot(norm) * tan.jhat /norm.jhat; u = u / (tan.ihat - (norm.ihat * tan.jhat / norm.jhat)); v = - (u1.dot(norm) + u * norm.ihat) / norm.jhat; at(j = y[0]) { /* And construct the new state vector */ localug[i,j].a = localug[i,j+1].a; localug[i,j].b = rho * u; localug[i,j].c = rho * v; localug[i,j].d = rho * (Cv * localtg[i,j] + 0.5 * (u*u + v*v)); localpg[i,j] = localpg[i,j+1]; } at(j = y[jmax-1]) { /* Top Wall Boundary Cells */ /* Checked numerically for default conditions. */ /* Construct normal and tangent vectors */ /* This part checked and works; it produces the correct vectors */ tan.ihat = xnode[i,j] - xnode[i-1,j]; tan.jhat = ynode[i,j] - ynode[i-1,j]; norm.ihat = ynode[i,j] - ynode[i-1,j]; norm.jhat = -(xnode[i,j] - xnode[i-1,j]); } scrap = tan.magnitude(); tan.ihat = tan.ihat / scrap; tan.jhat = tan.jhat / scrap; scrap = norm.magnitude(); norm.ihat = norm.ihat / scrap; norm.jhat = norm.jhat / scrap; at(j = y[jmax-1]) { /* Now set some state variables */ rho = localug[i,j].a; temp = localtg[i,j]; u1.ihat = localug[i,j].b / rho; u1.jhat = localug[i,j].c / rho; } u = u1.dot(tan) + u1.dot(norm) * tan.jhat /norm.jhat; u = u / (tan.ihat - (norm.ihat * tan.jhat / norm.jhat)); v = - (u1.dot(norm) + u * norm.ihat) / norm.jhat; at(j = y[jmax]) { /* And construct the new state vector */ localug[i,j].a = localug[i,j-1].a; localug[i,j].b = rho * u; localug[i,j].c = rho * v; localug[i,j].d = rho * (Cv * temp + 0.5 * (u*u + v*v)); localtg[i,j] = temp; localpg[i,j] = localpg[i,j-1]; } } Adlib.writeHalo(localtg, wlo, whi); Adlib.writeHalo(localug, wlo, whi); Adlib.writeHalo(localpg, wlo, whi); overall(j = y for 1 : jmax - 1) { at(i = x[0]){ /* Inlet Boundary Cells: unchecked */ /* Construct the normal vector; */ norm.ihat = ynode[i,j-1] - ynode[i,j]; norm.jhat = xnode[i,j] - xnode[i,j-1]; scrap = norm.magnitude(); norm.ihat = norm.ihat / scrap; norm.jhat = norm.jhat / scrap; theta = Math.acos((ynode[i,j-1] - ynode[i,j]) / Math.sqrt((xnode[i,j] - xnode[i,j-1])* (xnode[i,j] - xnode[i,j-1]) + (ynode[i,j-1] - ynode[i,j]) * (ynode[i,j-1] - ynode[i,j]))); u1.ihat = localug[i+1,j].b / localug[i+1,j].a; u1.jhat = localug[i+1,j].c / localug[i+1,j].a; uprime = u1.ihat * Math.cos(theta); c = Math.sqrt(gamma * rgas * localtg[i+1,j]); /* Supersonic inflow; works on the initial conditions */ if (uprime < -c) { /* Use far field conditions */ localug[i,j].a = rhoff; localug[i,j].b = rhoff * uff; localug[i,j].c = rhoff * vff; localug[i,j].d = rhoff * (Cv * tff + 0.5 * (uff*uff + vff*vff)); localtg[i,j] = tff; localpg[i,j] = pff; } /* Subsonic inflow */ /* This works on the initial conditions */ else if(uprime < 0.0) { /* Calculate Riemann invarients here */ jminus = u1.ihat - 2.0/(gamma-1.0) * c; s = Math.log(pff) - gamma * Math.log(rhoff); v = vff; u = (jplusff + jminus) / 2.0; scrap = (jplusff - u) * (gamma-1.0) * 0.5; localtg[i,j] = (1.0 / (gamma * rgas)) * scrap * scrap; localpg[i,j] = Math.exp(s) / Math.pow((rgas *localtg[i,j]), gamma); localpg[i,j] = Math.pow(localpg[i,j], 1.0 / (1.0 - gamma)); /* And now: construct the new state vector */ localug[i,j].a = localpg[i,j] / (rgas * localtg[i,j]); localug[i,j].b = localug[i,j].a * u; localug[i,j].c = localug[i,j].a * v; localug[i,j].d = localug[i,j].a * (Cv * tff + 0.5 * (u*u + v*v)); } /* Other options: */ /* We should throw an exception here. */ else { System.err.println("You have outflow at the inlet, "+ "which is not allowed."); } /* Outlet Boundary Cells */ /* Construct the normal vector; works. */ norm.ihat = ynode[i,j] - ynode[i,j-1]; norm.jhat = xnode[i,j-1] - xnode[i,j]; scrap = norm.magnitude(); norm.ihat = norm.ihat / scrap; norm.jhat = norm.jhat / scrap; scrap = xnode[i,j-1] - xnode[i,j]; scrap2 = ynode[i,j] - ynode[i,j-1]; theta = Math.acos((ynode[i,j] - ynode[i,j-1]) / Math.sqrt(scrap*scrap + scrap2*scrap2)); } at(i = x[imax-1]){ u1.ihat = localug[i,j].b / localug[i,j].a; u1.jhat = localug[i,j].c / localug[i,j].a; uprime = u1.ihat * Math.cos(theta); c = Math.sqrt(gamma * rgas * localtg[i,j]); } /* Supersonic outflow; works for defaults conditions. */ if (uprime > c){ at(i = x[imax-1]){ /* Use a backward difference 2nd order derivative approximation */ /* To set values at exit */ localug[i+1,j].a = 2.0 * localug[i,j].a - localug[i,j].a; localug[i+1,j].b = 2.0 * localug[i,j].b - localug[i-1,j].b; localug[i+1,j].c = 2.0 * localug[i,j].c - localug[i-1,j].c; localug[i+1,j].d = 2.0 * localug[i,j].d - localug[i-1,j].d; localpg[i+1,j] = 2.0 * localpg[i,j] - localpg[i-1,j]; localtg[i+1,j] = 2.0 * localtg[i,j] - localtg[i-1,j]; } } /* Subsonic Outflow; works for defaults conditions. */ else if (uprime < c && uprime > 0) { at(i = x[imax-1]){ jplus = u1.ihat + 2.0/(gamma - 1) * c; v = localug[i,j].c / localug[i,j].a; s = Math.log(localpg[i,j]) - gamma * Math.log(localug[i,j].a); u = (jplus + jminusff) / 2.0; scrap =(jplus - u)* (gamma-1.0) * 0.5; localtg[i+1,j] = (1.0 / (gamma * rgas)) * scrap * scrap; localpg[i+1,j] = Math.exp(s) / Math.pow((rgas * localtg[i+1,j]), gamma); localpg[i+1,j] = Math.pow(localpg[i+1,j], 1.0 / (1.0-gamma)); rho = localpg[i+1,j]/ (rgas * localtg[i+1,j]); /* And now, construct the new state vector. */ localug[i+1,j].a = rho; localug[i+1,j].b = rho * u; localug[i+1,j].c = rho * v; localug[i+1,j].d = rho * (Cv * localtg[i+1,j] + 0.5 * (u*u + v*v)); } } /* Other cases that shouldn't have to be used. */ else if (uprime < -c) { at(i = x[0]){ /* Supersonic inflow. Use far field conditions */ localug[i,j].a = rhoff; localug[i,j].b = rhoff * uff; localug[i,j].c = rhoff * vff; localug[i,j].d = rhoff * (Cv * tff + 0.5 * (uff*uff + vff*vff)); localtg[i,j] = tff; localpg[i,j] = pff; } } /* Subsonic inflow. This works on the initial conditions. */ else if(uprime < 0.0) { at(i = x[0]){ /* Debug: throw exception here? */ /* Calculate Riemann invarients here. */ jminus = u1.ihat - 2.0/(gamma-1.0) * c; s = Math.log(pff) - gamma * Math.log(rhoff); v = vff; u = (jplusff + jminus) / 2.0; scrap = (jplusff - u)* (gamma-1.0) * 0.5; localtg[i,j] = (1.0 / (gamma * rgas)) * scrap * scrap; localpg[i,j] = Math.exp(s) / Math.pow((rgas * localtg[i,j]), gamma); localpg[i,j] = Math.pow(localpg[i,j], 1.0 / (1.0 - gamma)); /* And now: construct the new state vector. */ localug[i,j].a = localpg[i,j] / (rgas * localtg[i,j]); localug[i,j].b = localug[i,j].a * u; localug[i,j].c = localug[i,j].a * v; localug[i,j].d = localug[i,j].a * (Cv * tff + 0.5 * (u*u + v*v)); } } /* Other Options */ /* Debug: throw exception here? */ else { System.err.println("You have inflow at the outlet, "+ "which is not allowed."); } } /* Do something with corners to avoid division by zero errors */ /* What you do shouldn't matter */ at(i = x[0]) at(j = y[0]) localug[i,j] = localug[i+1,j]; at(i = x[imax]) at(j = y[0]) localug[i,j] = localug[i,j+1]; at(i = x[0]) at(j = y[jmax]) localug[i,j] = localug[i+1,j]; at(i = x[imax]) at(j = y[jmax]) localug[i,j] = localug[i,j-1]; } void drawEdge(Graphics g, int width, int height) { /* Draws the outside border of the grid (as opposed to the frame) */ int gridheight, /* height of drawing region for grid in pixels */ gridwidth; /* Width of same */ int i, j; int x1, x2, y1, y2; int xoffset, yoffset; double scale; gridheight = height; /* Leave room at bottom for labels */ gridwidth = width; scale = Math.min(gridwidth/(xmax-xmin), gridheight/(ymax-ymin)); xoffset = (int) (xmin * scale); yoffset = (int) (ymin * scale); g.setColor(Color.white); for (i = 0; i < imax-1; ++i) { x1 = (int) (scale * xnodeE[i]) - xoffset; y1 = gridheight - (int) (scale * ynodeE[i]) - yoffset; x2 = (int) (scale * xnodeE[i+1]) - xoffset; y2 = gridheight - (int) (scale * ynodeE[i+1]) - yoffset; g.drawLine(x1, y1, x2, y2); x1 = (int) (scale * xnodeW[i]) - xoffset; y1 = gridheight - (int) (scale * ynodeW[i]) - yoffset; x2 = (int) (scale * xnodeW[i+1]) - xoffset; y2 = gridheight - (int) (scale * ynodeW[i+1]) - yoffset; g.drawLine(x1, y1, x2, y2); } for (j =0; j < jmax-1; ++j) { x1 = (int) (scale * xnodeS[j]) - xoffset; y1 = gridheight - (int) (scale * ynodeS[j]) - yoffset; x2 = (int) (scale * xnodeS[j+1]) - xoffset; y2 = gridheight - (int) (scale * ynodeS[j+1]) - yoffset; g.drawLine(x1, y1, x2, y2); x1 = (int) (scale * xnodeN[j]) - xoffset; y1 = gridheight - (int) (scale * ynodeN[j]) - yoffset; x2 = (int) (scale * xnodeN[j+1]) - xoffset; y2 = gridheight - (int) (scale * ynodeN[j+1]) - yoffset; g.drawLine(x1, y1, x2, y2); } } void drawTufts(Graphics g, int width, int height) { Color tuftscolor; double scale, theta; int gridwidth, gridheight; int i, j; int x1, x2, y1, y2; int xoffset, yoffset; tuftscolor = new Color(150,150,150); gridheight = height; /* Leave room for labels if necessary */ gridwidth = width; g.setColor(tuftscolor); scale = Math.min(gridwidth/(xmax-xmin), gridheight/(ymax-ymin)); xoffset = (int) (xmin * scale); yoffset = (int) (ymin * scale); on(p) { overall(i = x for 0 : imax - 2) overall(j = y for 0 : jmax - 2) { x1 = (int) (scale * xnode[i,j]) - xoffset; y1 = gridheight - (int) (scale * ynode[i,j]) - yoffset; if (Math.abs(tempug[i,j].b) > 0.00001) { theta = -Math.atan(tempug[i,j].c/tempug[i,j].b); } else { if (tempug[i,j].c > 0.0) theta = -Math.PI * 0.5; else theta = -Math.PI * 1.5; } x2 = x1 + (int) (6.0*Math.cos(theta)); y2 = y1 + (int) (6.0*Math.sin(theta)); g.drawLine(x1, y1, x2, y2); } } } void drawGrid(Graphics g, int width, int height) { /* Draws the grid; does not draw outside edges */ Color gridcolor; int i, j; /* Dummy Counters */ int gridheight, /* height of drawing region for grid in pixels */ gridwidth; /* Width of same */ int x1, y1, x2, y2; int xoffset, yoffset; double scale; gridheight = height; /* Leave room for labels if necessary */ gridwidth = width; gridcolor = new Color(100,100,100); g.setColor(gridcolor); scale = Math.min(gridwidth/(xmax-xmin), gridheight/(ymax-ymin)); xoffset = (int) (xmin * scale); yoffset = (int) (ymin * scale); on(p) { overall(i = x for 1 : imax - 1) overall(j = y for 1 : jmax - 1) { x1 = (int) (scale * xnode[i-1,j-1]) - xoffset; y1 = gridheight - (int) (scale * ynode[i-1,j-1]) - yoffset; x2 = (int) (scale * xnode[i,j-1]) - xoffset; y2 = gridheight - (int) (scale * ynode[i,j-1]) - yoffset; g.drawLine(x1, y1, x2, y2); x2 = (int) (scale * xnode[i-1,j]) - xoffset; y2 = gridheight - (int) (scale * ynode[i-1,j]) - yoffset; g.drawLine(x1, y1, x2, y2); } overall(i = x for 1 : imax - 1 ) { at(j = y[jmax-1]) { x1 = (int) (scale * xnode[i-1,j]) - xoffset; y1 = gridheight - (int) (scale * ynode[i-1,j]) - yoffset; x2 = (int) (scale * xnode[i,j]) - xoffset; y2 = gridheight - (int) (scale * ynode[i,j]) - yoffset; g.drawLine(x1, y1, x2, y2); } } at(i = x[imax-1]) { overall(j = y for 1 : jmax - 1) { x1 = (int) (scale * xnode[i,j-1]) - xoffset; y1 = gridheight - (int) (scale * ynode[i,j-1]) - yoffset; x2 = (int) (scale * xnode[i,j]) - xoffset; y2 = gridheight - (int) (scale * ynode[i,j]) - yoffset; g.drawLine(x1, y1, x2, y2); } } } } PlotData[] plotdata; synchronized void gatherLinePlotdata(int ndata){ double temp, vel, c; double u, v; double mach, pres; int j; plotdata = new PlotData[3]; for(int n = 0; n < 3; n++){ if (n == 0){ j = 1; }else if (n==1){ j = (jmax + 1) / 2; }else j = jmax; double[[-]] vals = new double[[x]] on p/y[j]; on(p){ at(jj = y[j]) overall(i = x for 1 : imax){ switch(ndata) { case 0: /* Density */ vals[i] = ug[i,jj].a; break; case 1: /* Pressure */ u = ug[i,jj].b/ug[i,jj].a; v = ug[i,jj].c/ug[i,jj].a; vals[i] = ug[i,jj].d/ug[i,jj].a - (0.5 * (u*u + v*v)); vals[i] = vals[i]/Cv; vals[i] = vals[i] * rgas * ug[i,jj].a; break; case 2: /* Temperature */ u = ug[i,jj].b/ug[i,jj].a; v = ug[i,jj].c/ug[i,jj].a; vals[i] = ug[i,jj].d/ug[i,jj].a - (0.5 * (u*u + v*v)); vals[i] = vals[i]/Cv; break; case 3: /* Mach Number */ u = ug[i,jj].b / ug[i,jj].a; v = ug[i,jj].c / ug[i,jj].a; vel = Math.sqrt(u*u + v*v); temp = ug[i,jj].d/ug[i,jj].a - (0.5 * (u*u + v*v)); temp = temp/Cv; c = Math.sqrt(gamma * rgas * temp); vals[i] = vel / c; break; case 4: /* Axial Velocity */ vals[i] = ug[i,jj].b / ug[i,jj].a; break; case 5: /* Normal Velocity */ vals[i] = ug[i,jj].c / ug[i,jj].a; break; case 6: /* Stagnation (Total) Pressure */ u = ug[i,jj].b / ug[i,jj].a; v = ug[i,jj].c / ug[i,jj].a; vel = Math.sqrt(u*u + v*v); temp = ug[i,jj].d/ug[i,jj].a - (0.5 * (u*u + v*v)); temp = temp/Cv; c = Math.sqrt(gamma * rgas * temp); mach = vel/c; mach = (1.0 + (gamma-1)*0.5*mach*mach); pres = ug[i,jj].d/ug[i,jj].a - (0.5 * (u*u + v*v)); pres = pres/Cv; pres = pres * rgas * ug[i,jj].a; vals[i] = pres * Math.pow(mach, gamma/(gamma-1.0)); break; default: break; } /* ends the switch */ } } plotdata[n] = new PlotData(Adlib.maxval(vals), Adlib.minval(vals), vals); Adlib.writeHalo(plotdata[n].values, wlo, whi); } } void drawLinePlot(Graphics g,int ndata, int width, int height) { int gridheight, gridwidth; int jmid; int plotheight, plotwidth; double scalex, scaley; int xoffset, yoffset; int n; int x1, y1, x2, y2; int xmargin = 55; int ymargin = 30; gridheight = height; gridwidth = width; plotheight = height-ymargin-20; /* Leave room at bottom for labels */ plotwidth = width-xmargin-20; jmid = (jmax + 1) / 2; double plotmax = plotdata[0].maxval > plotdata[1].maxval ? plotdata[0].maxval : plotdata[1].maxval; plotmax = plotmax > plotdata[2].maxval ? plotmax : plotdata[2].maxval; double plotmin = plotdata[0].minval < plotdata[1].minval ? plotdata[0].minval : plotdata[1].minval; plotmin = plotmin < plotdata[2].minval ? plotmin : plotdata[2].minval; /* Now draw it up! */ /* Start with Axes */ g.setColor(Color.white); if (plotmin > 0 && plotmax > 0) { plotmin = 0; plotmax *=1.20; } else if (plotmin < 0 && plotmax < 0) { plotmax = 0; plotmin *= 1.20; } /* Draw x axis */ int yval = gridheight + (int) ( (plotmin / (plotmax - plotmin)) * (double) plotheight)-ymargin; g.drawLine(xmargin, yval, xmargin+plotwidth, yval); /* Draw y axis */ g.drawLine(xmargin, gridheight-ymargin, xmargin, gridheight-ymargin-plotheight); scalex = (double) plotwidth/(xmax-xmin); scaley = (double) plotheight/(plotmax-plotmin); xoffset = (int) (xmin * scalex); yoffset = (int) (plotmin * scaley); if (ndata == 3) { g.setColor(new Color(128, 128, 128)); yval = gridheight -ymargin - (int)(1.0*scaley)+yoffset; g.drawLine(xmargin, yval, xmargin+plotwidth, yval); g.setColor(Color.white); g.drawString("1.0", 5, yval); } StringBuffer writeBuffer = new StringBuffer().append(plotmin); writeBuffer.setLength(4); g.drawString(writeBuffer.toString(), 5, gridheight-ymargin); g.drawString(String.valueOf(xmin), xmargin, gridheight-5); StringBuffer writeBuffer2 = new StringBuffer().append(plotmax); writeBuffer2.setLength(4); g.drawString(writeBuffer2.toString(), 5, gridheight-plotheight-ymargin+5); g.drawString(String.valueOf(xmax), plotwidth+xmargin-10, gridheight-5); x1 = gridwidth/2 - 10; g.drawString("Axial Position(x)", x1, gridheight-5); g.drawString("Red: Top Wall", xmargin+5, gridheight-(ymargin+25)); g.drawString("Blue: Centerline", xmargin+5, gridheight-(ymargin+15)); g.drawString("Green: Lower Wall", xmargin+5, gridheight-(ymargin+5)); g.setColor(Color.green); on(p){ overall(i = x for 1 : imax - 1) at(j = y[0]){ x1 = (int) (scalex * xnode[i-1,j]) - xoffset + xmargin; y1 = gridheight - ymargin - (int) (plotdata[0].values[i] * scaley)+ yoffset; x2 = (int) (scalex * xnode[i,j]) - xoffset + xmargin; y2 = gridheight - ymargin - (int)(plotdata[0].values[i+1]* scaley) + yoffset; g.drawLine(x1, y1, x2, y2); g.drawOval(x1-3, y1-3, 6,6); } g.setColor(Color.blue); overall(i = x for 1 : imax - 1) at(j = y[jmid]) { x1 = (int) (scalex * xnode[i-1,j]) - xoffset + xmargin; y1 = gridheight - ymargin - (int) (plotdata[1].values[i] * scaley) + yoffset; x2 = (int) (scalex * xnode[i,j]) - xoffset + xmargin; y2 = gridheight - ymargin - (int) (plotdata[1].values[i+1] * scaley) + yoffset; g.drawLine(x1, y1, x2, y2); g.drawOval(x1-3, y1-3, 6,6); } g.setColor(Color.red); overall(i = x for 1 : imax - 1) at(j = y[jmax-1]) { x1 = (int) (scalex * xnode[i-1,j]) - xoffset + xmargin; y1 = gridheight - ymargin - (int) (plotdata[2].values[i] * scaley) + yoffset; x2 = (int) (scalex * xnode[i,j]) - xoffset + xmargin; y2 = gridheight - ymargin - (int) (plotdata[2].values[i+1] * scaley) + yoffset; g.drawLine(x1, y1, x2, y2); g.drawOval(x1-3, y1-3, 6,6); } } } PlotData data; synchronized void gatherContourdata(int ndata){ HPutil.copy(tempug, ug); double temp, vel, c; double u, v; double mach, pres; double[[-,-]] vals = new double[[x,y]] on p; on(p) { overall(i = x for 1 : imax) overall(j = y for 1 : jmax){ switch(ndata){ case 0: /* Density */ vals[i,j] = ug[i,j].a; break; case 1: /* Pressure */ u = ug[i,j].b/ug[i,j].a; v = ug[i,j].c/ug[i,j].a; vals[i,j] = ug[i,j].d/ug[i,j].a - (0.5 * (u*u + v*v)); vals[i,j] = vals[i,j]/Cv; vals[i,j] = vals[i,j] * rgas * ug[i,j].a; break; case 2: /* Temperature */ u = ug[i,j].b/ug[i,j].a; v = ug[i,j].c/ug[i,j].a; vals[i,j] = ug[i,j].d/ug[i,j].a - (0.5 * (u*u + v*v)); vals[i,j] = vals[i,j]/Cv; break; case 3: /* Mach Number */ u = ug[i,j].b/ ug[i,j].a; v = ug[i,j].c / ug[i,j].a; vel = Math.sqrt(u*u + v*v); temp = ug[i,j].d/ug[i,j].a - (0.5 * (u*u + v*v)); temp = temp/Cv; c = Math.sqrt(gamma * rgas * temp); vals[i,j] = vel / c; break; case 4: /* Axial Velocity */ vals[i,j] = ug[i,j].b / ug[i,j].a; break; case 5: /* Normal Velocity */ vals[i,j] = ug[i,j].c / ug[i,j].a; break; case 6: /* Stagnation (Total) Pressure */ u = ug[i,j].b/ ug[i,j].a; v = ug[i,j].c / ug[i,j].a; vel = Math.sqrt(u*u + v*v); temp = ug[i,j].d/ug[i,j].a - (0.5 * (u*u + v*v)); temp = temp/Cv; c = Math.sqrt(gamma * rgas * temp); mach = vel/c; mach = (1.0 + (gamma-1)*0.5*mach*mach); pres = ug[i,j].d/ug[i,j].a - (0.5 * (u*u + v*v)); pres = pres/Cv; pres = pres * rgas * ug[i,j].a; vals[i,j] = pres * Math.pow(mach, gamma/(gamma-1.0)); break; default: /* Needs error checking here too?? */ break; } } } data = new PlotData(Adlib.maxval(vals), Adlib.minval(vals), vals); Adlib.writeHalo(data.data, wlo, whi); } void drawContour(Color ctable[], Graphics g, int ndata, int width, int height) { /* Draws contours outlining data set specified by ndata */ boolean drawcont=false; double c; /* Speed of sound */ double contval[]; double frac; double scale; double scrap; double uprime, vprime; double xnode1, xnode2; double ynode1, ynode2; int i1=0, i2=0, i3=0, i4=0; int j1=0, j2=0, j3=0, j4=0; int n; int gridheight, /* height of drawing region for grid in pixels */ gridwidth; /* Width of same */ int speciesindex; int xdraw1, xdraw2; int ydraw1, ydraw2; int xoffset, yoffset; /* Calculate the data you want from the state vector */ datamin = data.minval; datamax = data.maxval; /* Determine contour values */ contval = new double[ncontour]; scrap = datamax-datamin; for (n = 0; n < ncontour; ++n) contval[n] = datamin + ((double)n/(double)(ncontour-1))*scrap; /* Set screen values */ gridheight = height; /* Leave room at bottom for labels */ gridwidth = width; scale = Math.min(gridwidth/(xmax-xmin), gridheight/(ymax-ymin)); xoffset = (int) (xmin * scale); yoffset = (int) (ymin * scale); /* Go through cells looking for crossings */ on(p){ overall(i = x for 1 : imax - 1) overall(j = y for 1 : jmax - 1) { for (n = 0; n < ncontour; ++n) { /* First check the SW triangle */ drawcont = false; if ((contval[n] > data.data[i,j] && contval[n] < data.data[i+1,j]) || (contval[n] < data.data[i,j] && contval[n] > data.data[i+1,j])) { i1 = i`; j1 = j`; i2 = i`+1; j2 = j`; drawcont = true; } if ((contval[n] > data.data[i,j] && contval[n] < data.data[i,j+1]) || (contval[n] < data.data[i,j] && contval[n] > data.data[i,j+1])) if(drawcont) { i3=i`; j3=j`; i4 = i`; j4 = j`+1; } else { i1 = i`; j1 = j`; i2 = i`; j2 = j`+1; drawcont = true; } if ((contval[n] > data.data[i,j+1] && contval[n] < data.data[i+1,j]) || (contval[n] < data.data[i,j+1] && contval[n] > data.data[i+1,j])) { i3 = i`; j3 = j`+1; i4 = i`+1; j4 = j`; } if(drawcont) { /* Go ahead, draw the contour */ /* Sets color of contour */ i1 = i1 - i`; i2 = i2 - i`; i3 = i3 - i`; i4 = i4 - i`; j1 = j1 - j`; j2 = j2 - j`; j3 = j3 - j`; j4 = j4 - j`; g.setColor(ctable[n]); frac = Math.abs((contval[n]-data.data[i+i1,j+j1])/ (data.data[i+i1,j+j1] - data.data[i+i2,j+j2])); xnode1 = (xnode[i+(i2-1),j+(j2-1)]- xnode[i+(i1-1),j+(j1-1)]) * frac + xnode[i+(i1-1),j+(j1-1)]; xdraw1 = (int) (scale * xnode1) - xoffset; ynode1 = (ynode[i+(i2-1),j+(j2-1)] - ynode[i+(i1-1),j+(j1-1)]) * frac + ynode[i+(i1-1),j+(j1-1)]; ydraw1 = gridheight - (int) (scale * ynode1) - yoffset; frac = Math.abs((contval[n]-data.data[i+i3,j+j3])/ (data.data[i+i3,j+j3] - data.data[i+i4,j+j4])); xnode2 = (xnode[i+(i4-1),j+(j4-1)] - xnode[i+(i3-1),j+(j3-1)]) * frac + xnode[i+(i3-1),j+(j3-1)]; xdraw2 = (int) (scale * xnode2) - xoffset; ynode2 = (ynode[i+(i4-1),j+(j4-1)] - ynode[i+(i3-1),j+(j3-1)]) * frac + ynode[i+(i3-1),j+(j3-1)]; ydraw2 = gridheight - (int) (scale * ynode2) - yoffset; g.drawLine(xdraw1, ydraw1, xdraw2, ydraw2); } /* Then check the NE triangle */ drawcont = false; if ((contval[n] > data.data[i,j+1] && contval[n] < data.data[i+1,j+1]) || (contval[n] < data.data[i,j+1] && contval[n] > data.data[i+1,j+1])) { i1 = i`; j1 = j`+1; i2 = i`+1; j2 = j`+1; drawcont = true; } if ((contval[n] > data.data[i+1,j] && contval[n] < data.data[i+1,j+1]) || (contval[n] < data.data[i+1,j] && contval[n] > data.data[i+1,j+1])) if(drawcont) { i3=i`+1; j3=j`; i4 = i`+1; j4 = j`+1; } else { i1 = i`+1; j1 = j`; i2 = i`+1; j2 = j`+1; drawcont = true; } if ((contval[n] > data.data[i,j+1] && contval[n] < data.data[i+1,j]) || (contval[n] < data.data[i,j+1] && contval[n] > data.data[i+1,j])) { i3 = i`; j3 = j`+1; i4 = i`+1; j4 = j`; } if(drawcont) { /* Go ahead, draw the contour */ i1 = i1 - i`; i2 = i2 - i`; i3 = i3 - i`; i4 = i4 - i`; j1 = j1 - j`; j2 = j2 - j`; j3 = j3 - j`; j4 = j4 - j`; g.setColor(ctable[n]); frac = Math.abs((contval[n]-data.data[i+i1,j+j1])/ (data.data[i+i1,j+j1] - data.data[i+i2,j+j2])); xnode1 = (xnode[i+(i2-1),j+(j2-1)]- xnode[i+(i1-1),j+(j1-1)]) * frac + xnode[i+(i1-1),j+(j1-1)]; xdraw1 = (int) (scale * xnode1) - xoffset; ynode1 = (ynode[i+(i2-1),j+(j2-1)]- ynode[i+(i1-1),j+(j1-1)]) * frac + ynode[i+(i1-1),j+(j1-1)]; ydraw1 = gridheight - (int) (scale * ynode1) - yoffset; frac = Math.abs((contval[n]-data.data[i+i3,j+j3])/ (data.data[i+i3,j+j3] - data.data[i+i4,j+j4])); xnode2 = (xnode[i+(i4-1),j+(j4-1)]- xnode[i+(i3-1),j+(j3-1)]) * frac + xnode[i+(i3-1),j+(j3-1)]; xdraw2 = (int) (scale * xnode2) - xoffset; ynode2 = (ynode[i+(i4-1),j+(j4-1)] - ynode[i+(i3-1),j+(j3-1)]) * frac + ynode[i+(i3-1),j+(j3-1)]; ydraw2 = gridheight - (int) (scale * ynode2) - yoffset; g.drawLine(xdraw1, ydraw1, xdraw2, ydraw2); } } /* Ends "for (n = 0; n < ncoutour...)" */ } } } } class PlotData{ double maxval, minval; double [[-]] values; double [[-,-]] data; PlotData(double maxval, double minval, double[[-]] values){ this.maxval = maxval; this.minval = minval; this.values = values; } PlotData(double maxval, double minval, double[[-,-]] data){ this.maxval = maxval; this.minval = minval; this.data = data; } } class Statevector implements Serializable{ double a; /* Storage for Statevectors */ double b; double c; double d; Statevector() { a = 0.0; b = 0.0; c = 0.0; d = 0.0; } /* Most of these vector manipulation routines are not */ /* used in this program because I inlined them for speed. */ /* I leave them here because they may be useful in the future. */ public Statevector amvect(double m, Statevector that) { /* Adds statevectors multiplies the sum by scalar m */ Statevector answer = new Statevector(); answer.a = m * (this.a + that.a); answer.b = m * (this.b + that.b); answer.c = m * (this.c + that.c); answer.d = m * (this.d + that.d); return answer; } public Statevector avect(Statevector that) { Statevector answer = new Statevector(); /* Adds two statevectors */ answer.a = this.a + that.a; answer.b = this.b + that.b; answer.c = this.c + that.c; answer.d = this.d + that.d; return answer; } public Statevector mvect(double m) { Statevector answer = new Statevector(); /* Multiplies statevector scalar m */ answer.a = m * this.a; answer.b = m * this.b; answer.c = m * this.c; answer.d = m * this.d; return answer; } public Statevector svect(Statevector that) { Statevector answer = new Statevector(); /* Subtracts vector that from this */ answer.a = this.a - that.a; answer.b = this.b - that.b; answer.c = this.c - that.c; answer.d = this.d - that.d; return answer; } public Statevector smvect(double m, Statevector that) { Statevector answer = new Statevector(); /* Subtracts statevector that from this and multiplies the */ /* result by scalar m */ answer.a = m * (this.a - that.a); answer.b = m * (this.b - that.b); answer.c = m * (this.c - that.c); answer.d = m * (this.d - that.d); return answer; } } class Vector2 { double ihat; /* Storage for 2-D vector */ double jhat; Vector2() { ihat = 0.0; jhat = 0.0; } public double magnitude() { double mag; mag = Math.sqrt(this.ihat*this.ihat + this.jhat * this.jhat); return mag; } public double dot(Vector2 that) { /* Calculates dot product of two 2-d vector */ double answer; answer = this.ihat * that.ihat + this.jhat * that.jhat; return answer; } } class DataPanel extends HPPanel implements HPspmd{ boolean gridSwitch = false; boolean tuftsSwitch = false; Color ctable[]; datafield data2D; public double machff; public double fourthOrderDamping; public double secondOrderDamping; public int nplot; public int ntime; public int nview; volatile int itercount; final int iterationsPerCycle = 1; /* Number of iterations per redraw */ final int keyWidth=100; /* width of color key in pixels */ /* For the screen buffering */ private Image offScreenImage; private Graphics offScreenGraphics; private java.awt.Dimension offScreenSize; DataPanel(InputStream instream, Procs control) {/* Read in the input file */ itercount = 0; try { data2D = new datafield(instream, control); } catch (IOException e) { System.err.println("Input/output error: probably a bad "+ "filename or bad connection." + e.getMessage()); } machff = data2D.machfftarget; fourthOrderDamping = 1.0; secondOrderDamping = 1.0; /* Initialize Color Table */ ctable = new Color[16]; ctable[0] = new Color(0,0,255); ctable[1] = new Color(0,64,255); ctable[2] = new Color(0,128,255); ctable[3] = new Color(0,192,255); ctable[4] = new Color(0,255,255); ctable[5] = new Color(0,255,192); ctable[6] = new Color(0,255,128); ctable[7] = new Color(0,255,64); ctable[8] = new Color(0,255,0); ctable[9] = new Color(32,255,0); ctable[10] = new Color(96,255,0); ctable[11] = new Color(160,255,0); ctable[12] = new Color(255,255,0); ctable[13] = new Color(255,160,0); ctable[14] = new Color(255,96,0); ctable[15] = new Color(255,0,0); } public void broadcastInitialData(){ data2D.broadcastInitialData(); } public void doIteration() { int n; itercount++; for (n = 0; n < iterationsPerCycle; ++n) { data2D.machfftarget = machff; data2D.secondOrderDamping = secondOrderDamping; data2D.fourthOrderDamping = fourthOrderDamping; data2D.ntime = ntime; data2D.doIteration(); } if (nplot == 0) { data2D.gatherContourdata(nview); } else data2D.gatherLinePlotdata(nview); repaint(); } public void setGridSwitch(boolean nchoice) { gridSwitch = nchoice; } public void setTuftsSwitch(boolean nchoice) { tuftsSwitch = nchoice; } public void hpPaint(Graphics g) { if (itercount == 0) return; g.setColor(Color.black); if (nplot == 0) { g.fillRect(0, 0, getSize().width-keyWidth, getSize().height - 1); if (gridSwitch) data2D.drawGrid(g, getSize().width-keyWidth, getSize().height-1); if (tuftsSwitch) data2D.drawTufts(g, getSize().width-keyWidth, getSize().height-1); data2D.drawContour(ctable, g, nview, getSize().width-keyWidth, getSize().height-1); data2D.drawEdge(g, getSize().width-keyWidth, getSize().height-1); data2D.drawKey(ctable, g, keyWidth, getSize().height-1, getSize().width-1, getSize().height-1); } else { g.fillRect(0, 0, getSize().width-1, getSize().height - 1); data2D.drawLinePlot(g, nview, getSize().width-1, getSize().height-1); } } boolean linePlotcalled = false; public final synchronized void hpUpdate (Graphics theG) { /* This section is based on a borrowed (Matt Gray - MIT) */ /* NoFlickerApplet. (impliments Double Buffering) */ java.awt.Dimension d = getSize(); if((offScreenImage == null) || (d.width != offScreenSize.width) || (d.height != offScreenSize.height)) { offScreenImage = createImage(d.width, d.height); offScreenSize = d; offScreenGraphics = offScreenImage.getGraphics(); offScreenGraphics.setFont(getFont()); } offScreenGraphics.setColor(Color.black); if (nplot == 0) { if (linePlotcalled) { offScreenGraphics.fillRect(0, 0, getSize().width-1, getSize().height - 1); linePlotcalled = false; } else offScreenGraphics.fillRect(0,0,getSize().width-keyWidth, getSize().height - 1); if (gridSwitch) data2D.drawGrid(offScreenGraphics, getSize().width-keyWidth, getSize().height-1); if (tuftsSwitch) data2D.drawTufts(offScreenGraphics, getSize().width-keyWidth, getSize().height-1); data2D.drawContour(ctable, offScreenGraphics, nview, getSize().width-keyWidth, getSize().height-1); data2D.drawEdge(offScreenGraphics, getSize().width-keyWidth, getSize().height-1); data2D.drawKey(ctable, offScreenGraphics, keyWidth, getSize().height-1, getSize().width-1, getSize().height-1); } else { linePlotcalled = true; offScreenGraphics.fillRect(0, 0, getSize().width-1, getSize().height - 1); data2D.drawLinePlot(offScreenGraphics, nview, getSize().width-1, getSize().height-1); } theG.drawImage(offScreenImage, 0, 0, null); } } class ParameterPanel extends Panel { final int sliderMin = 0; final int sliderMax = 100; double sliderValMin; double sliderValue; double sliderValMax; Scrollbar slider; TextField textField; double getValue() { double f; try { f = Double.valueOf(textField.getText()).doubleValue(); } catch (java.lang.NumberFormatException e) { f = sliderValMin; } return f; } void makeSlidePanel(String parameterTitle) { /* This lays out components, sets initial appearances */ GridBagConstraints c = new GridBagConstraints(); GridBagLayout gridbag = new GridBagLayout(); setLayout(gridbag); c.fill = GridBagConstraints.HORIZONTAL; Label label = new Label (parameterTitle, Label.LEFT); c.weightx = 0.0; gridbag.setConstraints(label, c); add(label); textField = new TextField(String.valueOf(sliderValue), 6); c.weightx = 1.0; c.gridwidth = GridBagConstraints.REMAINDER; gridbag.setConstraints(textField, c); add(textField); int sliderIntValue = (int) ((sliderValue-sliderValMin)/(sliderValMax-sliderValMin) * (double) (sliderMax-sliderMin)) + sliderMin; slider = new Scrollbar(Scrollbar.HORIZONTAL, sliderIntValue, 10, sliderMin, sliderMax); c.weightx = 0.0; c.gridheight = 1; c.gridwidth = GridBagConstraints.REMAINDER; gridbag.setConstraints(slider, c); add(slider); } void setValue(int n) { setSliderValue(n); textField.setText(String.valueOf(sliderValue)); } void setSliderValue(double f) { sliderValue = f; int sliderIntValue; sliderIntValue = (int) ((f-sliderValMin) / (sliderValMax-sliderValMin) * (double) (sliderMax-sliderMin)) + sliderMin; if (sliderIntValue > sliderMax) { sliderIntValue = sliderMax; sliderValue = sliderValMax; } if (sliderIntValue < sliderMin) { sliderIntValue = sliderMin; sliderValue = sliderValMin; } slider.setValue(sliderIntValue); textField.setText(String.valueOf(sliderValue)); } public boolean handleEvent(Event e) { if (e.target instanceof Scrollbar) { sliderValue = ((double) slider.getValue()/ (double) (sliderMax-sliderMin)) * (sliderValMax-sliderValMin) + sliderValMin; textField.setText(String.valueOf(sliderValue)); return true; } else if (e.target instanceof TextField) { setSliderValue(getValue()); return true; } return false; } } class MachffPanel extends ParameterPanel { /* Ok, here's where we handle events for this scrollbar */ MachffPanel(String parameterTitle, double min, double max, double initialMachff) { super(); sliderValMin = min; sliderValMax = max; sliderValue = initialMachff; makeSlidePanel(parameterTitle); } public boolean handleEvent(Event e) { double temp; /* Overrides ParameterPanel's default behavior */ if (e.target instanceof Scrollbar) { sliderValue = ((double) slider.getValue()/ (double) (sliderMax-sliderMin)) * (sliderValMax-sliderValMin) + sliderValMin; if (sliderValue > 0.9 && sliderValue < 1.3) sliderValue = 1.3; textField.setText(String.valueOf(sliderValue)); return true; } else if (e.target instanceof TextField) { temp = getValue(); if (temp > 0.9 && temp < 1.3) { temp = 1.3; } setSliderValue(temp); return true; } return false; } } class DampingPanel extends ParameterPanel { /* Ok, here's where we handle events for this scrollbar */ DampingPanel(String parameterTitle, double min, double max, double initialDamping) { super(); sliderValMin = min; sliderValMax = max; sliderValue = initialDamping; makeSlidePanel(parameterTitle); } }