program main !*********************************************************************** ! !! MAIN is the main program for DISPLAY3. ! ! Discussion: ! ! DISPLAY3 reads graphics data from FLOW programs and makes plots. ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! "Primitive" parameters, that depend on nothing else. ! integer, parameter :: maxbou = 500 integer, parameter :: maxnpe = 6 integer, parameter :: maxnx = 81 integer, parameter :: maxny = 51 integer, parameter :: maxobj = 37 integer, parameter :: maxpar = 15 ! ! Parameters that depend on primitive parameters. ! integer, parameter :: maxelm = 2*(maxnx-1)*(maxny-1) integer, parameter :: maxnp = (2*maxnx-1)*(2*maxny-1) integer, parameter :: maxsen = 2*maxpar ! character ( len = 10 ) arrow real bval character ( len = 6 ) chrint character ( len = 80 ) command real cp(maxnp) real delx real dely character ( len = 10 ) dev real dudxn(maxnp) real dudyn(maxnp) real dvdxn(maxnp) real dvdyn(maxnp) logical echo logical eflag(maxelm) logical eflagu(maxelm) character ( len = 2 ) eqn(3,maxnp) real etaref(maxnpe) character ( len = 80 ) fildat character ( len = 80 ) filgrf character ( len = 80 ) filelm character ( len = 80 ) filinp character ( len = 20 ) filtyp real grace real gval integer i integer icmax integer icmin integer icolor(maxobj) integer icomp integer idata integer ierror integer ifile logical inside integer iplot character isay integer iset integer isotri(maxelm) integer itable integer itemp integer iwork1(maxelm) integer iwork2(maxnp) integer iwrite integer j integer jbound(5,maxbou) integer jcmax integer jcmin integer jfile integer jtemp integer k character ( len = 30 ) labelx character ( len = 30 ) labely logical lbar integer lenc integer lens integer lent integer line(maxobj) logical lppro logical lptpro logical lupro logical lutpro logical lvpro logical lvtpro integer nbound integer ncon integer nelem logical nflag(maxnp) logical nflag0(maxnp) integer node(maxnpe,maxelm) logical none integer np integer npar integer npe integer nprof(2*maxny-1) integer nsen integer nset integer nx integer nxskip integer ny integer nyskip character ( len = 30 ) object(maxobj) logical ovrlay real p(maxnp,0:maxsen) real para(maxpar) real ptar(maxnp) real rho(maxnp) real rmach(maxnp) real rval real s(maxnp) logical s_eqi real s2(maxnp) real scalee real scalen real scalev logical show(maxobj) real smax real smin real srange character ( len = 80 ) string real t(maxnp) real temp character ( len = 40 ) title character ( len = 40 ) title2 real u(maxnp,0:maxsen) real utar(maxnp) real v(maxnp,0:maxsen) real vtar(maxnp) real x1max real x1maxc real x1min real x1minc real x2max real x2min real x4max real x4min real xc(maxnp) real xmax real xmin real xprof real xsiref(maxnpe) real xsmax real xsmin real xtmax real xtmin real y1max real y1maxc real y1min real y1minc real y2max real y2min real y4max real y4min real yc(maxnp) real ymax real ymin real ysmax real ysmin real ytmax real ytmin ! call timestamp ( ) ! ! Greetings! ! call hello ( maxbou, maxelm, maxnp, maxnpe, maxnx, maxny, maxobj, & maxpar, maxsen ) ! ! Set initial values. ! call init ( arrow, cp, delx, dely, dev, dudxn, dudyn, dvdxn, dvdyn, echo, & eflag, eflagu, eqn, etaref, fildat, filgrf,filinp,filtyp,grace,icmax,icmin, & icolor,icomp,idata,ifile,iplot,iset,itable,iwrite,jcmax, & jcmin,labelx,labely,lbar,line,lppro,lptpro,lupro, & lutpro,lvpro,lvtpro,maxelm,maxnp,maxnpe,maxny,maxobj, & maxpar,maxsen,nbound,ncon,nelem,nflag,nflag0,node, & np,npe,nprof,nsen,nset,nxskip,ny,nyskip,object,ovrlay,p,rho,rmach,scalee, & scalen,scalev,show,smax,smin,title,title2,u,v,x1max,x1min, & x2max,x2min,x4max,x4min,xc,xprof,xsiref,xsmax,xsmin,y1max, & y1min,y2max,y2min,y4max,y4min,yc,ysmax,ysmin) ! ! Open the file to contain a copy of the user input commands. ! open ( unit = 17, file = filinp, status = 'replace', & access = 'sequential', & form = 'formatted' ) ! ! Get command from user ! 10 continue write(*,*)' ' write(*,*)'? ("H" for help)' read(*,'(a)',end = 80,err=80)command call s_blank_delete ( command ) write ( 17, '(a)' ) trim ( command ) if ( echo ) then write ( *, '(a)' ) trim ( command ) end if if ( command == ' ' ) then go to 10 end if ! ! A: advance, read next set of data from file. ! if ( s_eqi ( command,'a') ) then if ( s_eqi ( filtyp, 'FLOW' ) ) then call advanc ( delx,dely,echo,eflag,eqn,grace,ifile,iset,isotri, & jbound,maxbou,maxelm,maxnp,maxnpe,maxny,maxpar,maxsen, & nbound,nelem,nflag,nflag0,node,np,npar,npe,nprof,nsen, & nset,nx,ny,p,para,ptar,rho,srange,u,utar,v, & vtar,x1max,x1min,x2max,x2min,xc,xmax, & xmin,xsmax,xsmin,xtmax,xtmin,y1max,y1min,y2max,y2min,yc, & ymax,ymin,ysmax,ysmin,ytmax,ytmin) end if idata = 1 ! ! ARROWS = : Choose solid or line. ! else if ( s_eqi ( command(1:5),'arrow') ) then if ( s_eqi ( command(1:6),'arrow=') ) then arrow = command(7:) else if ( s_eqi ( command(1:7),'arrows=') ) then arrow = command(8:) else write(*,*)' ' write(*,*)'DISPLAY - Input request:' write(*,*)' Choose HOLLOW, LINE, SOLID for arrows.' read(*,'(a)')arrow write(17,'(a)')arrow if ( echo ) then write(*,'(a)')arrow end if end if write(*,*)' The arrow option is now ARROW = ' // trim ( arrow ) ! ! B: show boundary. ! else if ( s_eqi ( command, 'B' ) ) then show(1) = .not. show(1) if ( show(1) ) then write(*,*)'The boundary will be shown.' else write(*,*)'The boundary will NOT be shown.' end if ! ! BACK: show the background ! else if ( s_eqi ( command(1:4), 'BACK' ) ) then show(21) = .not.show(21) if ( show(21) ) then write(*,*)'The background will be shown.' else write(*,*)'The background will NOT be shown.' end if ! ! BAR: Switch display of color bar. ! else if ( s_eqi ( command, 'BAR' ) ) then lbar = .not.lbar if ( lbar ) then write(*,*)'The color bar will be shown.' else write(*,*)'The color bar will NOT be shown.' end if ! ! BH: bottom half. ! else if ( s_eqi ( command,'BH') ) then ysmax = ysmin+0.5*(ysmax-ysmin) call pltbox ( grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! BL: bottom left quarter. ! else if ( s_eqi ( command,'BL') ) then xsmax = xsmin+0.5*(xsmax-xsmin) ysmax = ysmin+0.5*(ysmax-ysmin) call pltbox ( grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! BC: bottom center quarter. ! else if ( s_eqi ( command,'BC') ) then temp = 0.25*(xsmax-xsmin) xsmin = xsmin+temp xsmax = xsmax-temp ysmax = ysmin+0.5*(ysmax-ysmin) call pltbox ( grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! BR: bottom right quarter. ! else if ( s_eqi ( command,'br') ) then xsmin = xsmin+0.5*(xsmax-xsmin) ysmax = ysmin+0.5*(ysmax-ysmin) call pltbox ( grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! C: choose colors. ! else if ( s_eqi ( command,'c') ) then write(*,*)' ' write(*,*)' Number Color Name' write(*,*)' ' do i = 1,maxobj write(*,'(1x,i2,2x,i3,2x,a)')i,icolor(i),object(i) end do write(*,*)' ' write(*,*)'Enter an object number, and a color number.' read(*,*,end = 90,err=90)itemp,jtemp write(17,*)itemp,jtemp if ( echo ) then write(*,*)itemp,jtemp end if if ( 1 <= itemp.and.itemp <= maxobj ) then icolor(itemp) = jtemp else write(*,*)'Your object number was out of bounds.' end if ! ! CC: choose color contour labels ! ! For some strange reason, in order to make the color table ! active, we have to call NEWFRM! ! else if ( s_eqi ( command(1:2),'cc') ) then if ( dev == ' ' ) then write(*,*)' ' write(*,*)'Please use the DEV command first!' go to 10 end if if ( s_eqi ( command(1:3),'cc=') ) then read(command(4:),*,end = 90,err=90)itable else write(*,*)' ' write(*,*)'Built in color tables include:' write(*,*)' ' write(*,*)'1 low black to high white.' write(*,*)'2 low blue to high yellow.' write(*,*)'3 low red, high blue, with bands between.' write(*,*)'4 low red, yellow, green, blue, high white.' write(*,*)'5 low white, blue, green, yellow, high red' write(*,*)'6 low blue to high red.' write(*,*)'7 linear table between 2 user colors.' write(*,*)'8 linear table between N user colors.' write(*,*)'9 low white to high black.' write(*,*)' ' write(*,*)'Enter a color table index between 1 and 9,' write(*,*)'or 0 to enter a color table from a file.' read(*,*,end = 90,err=90)itable write(17,*)itable if ( echo ) then write(*,*)itable end if end if call get_tab ( dev,echo,filgrf,grace,icmax,icmin,ierror,iplot,itable,ovrlay ) if ( itable == 1.or.itable == 9 ) then jcmax = 200 jcmin = 32 else jcmax = 255 jcmin = 2 end if write(*,*)' ' write(*,*)'Lowest color used will be JCMIN = ',jcmin write(*,*)'Highest color used will be JCMAX = ',jcmax ! ! CH: center half. ! else if ( s_eqi ( command,'ch') ) then temp = 0.25*(xsmax-xsmin) xsmin = xsmin+temp xsmax = xsmax-temp call pltbox ( grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! COLOR ! else if ( s_eqi ( command,'color') ) then write(*,*)'Enter the color index between 0 and 255' read(*,*,end = 90,err=90)i write(17,*)i if ( echo ) then write(*,*)i end if write(*,*)'Enter(R,G,B)' write(*,*)'Note: (0,0,0) is black, (1,1,1) is white!' read(*,*,end = 90,err=90)rval,gval,bval write(17,*)rval,gval,bval if ( echo ) then write(*,*)rval,gval,bval end if call setclr(i,bval,gval,rval) ! ! CP: draw CP contour lines. ! else if ( s_eqi ( command,'cp') ) then show(29) = .not.show(29) if ( show(29) ) then write(*,*)'CP contours will be plotted.' else write(*,*)'CP contours will NOT be plotted.' end if ! ! CPC: draw CP colors. ! else if ( s_eqi ( command,'cpc') ) then show(32) = .not.show(32) if ( show(32) ) then write(*,*)'CP colors will be plotted.' else write(*,*)'CP colors will NOT be plotted.' end if ! ! CTAB ! else if ( s_eqi ( command,'ctab') ) then call preplt ( dev,echo,filgrf,icmax,icmin,iplot,itable,ovrlay ) call cbox ( grace ) ! ! If the user hadn't already graphed anything, X1MIN = X1MAX=0 is ! possible. Set to a default. ! if ( x1min == x1max ) then x1min = 0.0 x1max = 1.0 else if ( y1min == y1max ) then y1min = 0.0 y1max = 1.0 end if call setwcd ( x1min,y1min,x1max,y1max,ierror ) call buzz ( dev, x1minc, x1maxc, y1minc, y1maxc ) ! ! DAT = : specify the data file to be read. ! else if ( s_eqi ( command(1:3),'DAT') ) then filtyp = 'flow' if ( s_eqi ( command(1:4),'dat=') ) then fildat = command(5:) else if ( ifile > 0 ) then close(unit = 2) write(*,*)'OpnFil is closing the old file '// trim ( fildat ) end if write(*,*)' ' write(*,*)'Enter the name of the new input data file:' read(*,'(a)',err = 90,end=90) fildat write(17,'(a)') trim ( fildat ) if ( echo ) then write(*,'(a)') trim ( fildat ) end if end if call opnfil ( fildat, ifile, iset ) if ( ifile == -2)go to 110 if ( ifile == -1)go to 100 iset = 0 20 continue call plot_file_read ( eqn,ierror,ifile,iset,isotri,maxelm,maxnp, & maxnpe,maxny,maxpar,maxsen,nelem,node,np,npar,npe,nprof, & nsen,nx,ny,p,para,ptar,rho,u,utar,v,vtar,xc,xprof,yc) if ( ierror == 0 ) then iset = iset+1 go to 20 end if nset = iset string = 'This file contains '//chrint(nset)//' data sets.' call s_blanks_delete ( string) write ( *, '(a)') trim ( string ) xsmin = 0.0 xsmax = 0.0 ysmin = 0.0 ysmax = 0.0 ! ! DB: debug output. ! else if ( s_eqi ( command,'db') ) then call node_data_print ( maxnp,maxsen,nflag,np,p,u,v,xc,yc ) ! ! 'DEV = ' Choose the graphics device. ! else if ( s_eqi ( command(1:3),'dev') ) then if ( dev /= ' ' ) then write(*,*)' ' write(*,*)'Display - Error!' write(*,*)' You have already chosen device '// trim ( dev ) write(*,*)' You may not change your mind!' go to 10 end if if ( s_eqi ( command(1:4),'dev=') ) then dev = trim ( command(5:) ) else write ( *, * ) ' ' write ( *, * ) 'Enter the graphics device desired.' write ( *, * ) ' ' write ( *, * ) ' CGMB output to a CGM binary file.' write ( *, * ) ' PS output to a PostScript file.' write ( *, * ) ' XWS output to an X window screen.' read(*,'(a)',end = 90, err=90 ) dev write ( 17, '(a)' ) dev if ( echo ) then write ( *, '(a)' ) dev end if end if if ( s_eqi ( dev(1:3),'cgm') ) then dev = 'cgmb' write(*,*)'Output will be to a CGM binary file "display.cgm".' else if ( s_eqi ( dev,'ps') ) then write(*,*)'Output will be to a PostScript file "display.ps".' else if ( s_eqi ( dev,'xws') ) then write(*,*)'Output will be to an X window screen.' else write(*,*)'Your device '// trim ( dev ) //' was not recognized!' dev = ' ' end if ! ! DJE = : specify the JEFF element file to be read. ! else if ( s_eqi ( command(1:3),'dje') ) then filtyp = 'jeff' if ( s_eqi ( command(1:4),'dje=') ) then filelm = command(5:) else write(*,*)'Enter the name of the element data file.' read(*,'(a)')filelm write(17,'(a)')filelm if ( echo ) then write(*,'(a)')filelm end if end if call rdelj(filelm,ifile,maxelm,maxnp,maxnpe,nelem,node,np,npe) xsmin = 0.0 xsmax = 0.0 ysmin = 0.0 ysmax = 0.0 ! ! DJN = : specify the JEFF node file to be read. ! else if ( s_eqi ( command(1:3),'djn') ) then if ( s_eqi ( command(1:4),'djn=' ) ) then fildat = command(5:) else write(*,*)'Enter the name of the node data file.' read(*,'(a)')fildat write(17,'(a)') trim ( fildat ) if ( echo ) then write(*,'(a)') trim ( fildat ) end if end if call adjeff(delx,dely,eflag,eqn,fildat,grace,ifile,isotri, & jbound,jfile,maxbou,maxelm,maxnp,maxnpe,nbound,nelem,nflag,nflag0, & node,np,npe,nset,nx,ny,p,rho,srange,u,v,x1max,x1min,x2max,x2min,xc, & xmax,xmin,xsmax,xsmin,xtmax,xtmin,y1max,y1min, & y2max,y2min,yc,ymax,ymin,ysmax,ysmin,ytmax,ytmin) ! ! DTEC = : specify the TECPLOT data file to be read. ! else if ( s_eqi ( command(1:4),'dtec') ) then filtyp = 'tecplot' if ( s_eqi ( command(1:5),'dtec=') ) then filelm = command(6:) else write(*,*)'Enter the name of the TECPLOT data file.' read(*,'(a)')filelm write(17,'(a)')filelm if ( echo ) then write(*,'(a)')filelm end if end if call rdtec(cp,filelm,ifile,maxelm,maxnp,maxnpe,nelem,node, & np,npe,nx,ny,rho,rmach,u,v,xc,yc) xsmin = 0.0 xsmax = 0.0 ysmin = 0.0 ysmax = 0.0 nset = 1 ifile = 2 ! ! We assume all triangles are nonisoperimetric. ! do i = 1,nelem isotri(i) = 0 end do ! ! Initially, all elements will be visible. ! do i = 1,nelem eflag(i) = .true. end do do i = 1,np nflag(i) = .true. nflag0(i) = .true. end do ! ! We give default values to the equations. ! do i = 1,np eqn(1,i) = 'U' eqn(2,i) = 'V' eqn(3,i) = 'P' end do ! ! Check element orientation. ! call element_check(maxnpe,nelem,node,np,npe,xc,yc) ! ! Determine the region size. ! call region_size ( delx,dely,grace,nelem,nflag,np,srange, & x1max,x1min,x2max,x2min,xc,xmax,xmin,xsmax,xsmin,xtmax,xtmin,y1max, & y1min,y2max,y2min,yc,ymax,ymin,ysmax,ysmin,ytmax,ytmin) ! ! Set the location of boundary edges. ! call boundary_set(eqn,jbound,maxbou,maxnpe,nbound,nelem,node,np,npe) ! ! E: show elements. ! else if ( s_eqi ( command,'e').or.s_eqi ( command,'elements') ) then show(2) = .not.show(2) if ( show(2) ) then write(*,*)'Elements will be shown.' else write(*,*)'Elements will NOT be shown.' end if ! ! EC: show element colors. ! else if ( s_eqi ( command,'ec') ) then show(20) = .not.show(20) if ( show(20) ) then write(*,*)'Element colors will be shown.' else write(*,*)'Element colors will NOT be shown.' end if ! ! EN: Show element numbers ! else if ( s_eqi ( command,'en') ) then show(28) = .not.show(28) if ( show(28) ) then write(*,*)'Element numbers will be shown.' else write(*,*)'Element numbers will NOT be shown.' end if ! ! EX: Define example problem. ! else if ( s_eqi ( command(1:2),'ex') ) then if ( ifile > 0 ) then close(unit = 2) write(*,*)'OpnFil is closing the old file '// trim ( fildat ) end if fildat = ' ' call exdat(eqn,isotri,maxelm,maxnp,maxnpe,maxny,maxpar, & maxsen,nelem,node,np,npar,nprof, & nx,ny,p,para,ptar,u,utar,v,vtar,xc,xprof,yc) do i = 1,np nflag(i) = .true. nflag0(i) = .true. end do do i = 1,nelem eflag(i) = .true. end do ! ! Determine the region size. ! call region_size ( delx,dely,grace,nelem,nflag,np,srange, & x1max,x1min,x2max,x2min,xc,xmax,xmin,xsmax,xsmin,xtmax,xtmin,y1max, & y1min,y2max,y2min,yc,ymax,ymin,ysmax,ysmin,ytmax,ytmin) call boundary_set(eqn,jbound,maxbou,maxnpe,nbound,nelem,node,np,npe) idata = 1 nset = 1 string = 'This example includes '//chrint(nset)//' data sets.' call s_blanks_delete ( string) write(*,'(1x,a)') trim ( string ) ! ! ECHO: switch the echo option ! else if ( s_eqi ( command,'echo') ) then echo = .not.echo if ( echo ) then write(*,'(a)') trim ( command ) write(*,*)'User input will be echoed.' else write(*,*)'User input will NOT be echoed.' end if ! ! F: switch the frame option ! else if ( s_eqi ( command,'f') ) then show(3) = .not.show(3) if ( show(3) ) then write(*,*)'The frame will be shown.' else write(*,*)'The frame will NOT be shown.' end if ! ! FILE = : set the name of the graphics output file. ! else if ( s_eqi ( command(1:4),'file') ) then if ( iplot > 0 ) then write(*,*)' ' write(*,*)'Display - Warning:' write(*,*)' It is too late to specify a plot file name.' else if ( s_eqi ( command(1:5),'file=') ) then filgrf = command(6:) else write(*,*)' ' write(*,*)'Enter the plot file name.' read(*,'(a)',err = 90,end=90)filgrf write(17,'(a)') trim ( filgrf ) if ( echo ) then write(*,'(a)') trim ( filgrf ) end if end if write(*,*)' ' write(*,*)'DISPLAY - Note:' write(*,*)' The plot file will be named '//filgrf ! ! FILTYP = : set the type of the input file. ! else if ( s_eqi ( command(1:6),'filtyp') ) then if ( s_eqi ( command(1:7),'filtyp=') ) then filtyp = command(8:) else write(*,*)' ' write(*,*)'Enter the file type (FLOW, JEFF or TECPLOT).' read(*,'(a)',err = 90,end=90)filtyp write(17,'(a)')filtyp if ( echo ) then write(*,'(a)')filtyp end if end if write(*,*)' ' write(*,*)'DISPLAY - Note:' write(*,*)' The input file type is set to '//filtyp ! ! FRAME: show the frame ! else if ( s_eqi ( command,'frame') ) then show(3) = .true. write(*,*)'DISPLAY - Note:' write(*,*)' A frame will be shown around the picture.' ! ! FULL: show full picture. ! else if ( s_eqi ( command,'full') ) then xsmax = xtmax xsmin = xtmin ysmax = ytmax ysmin = ytmin call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! G: create current graph. ! else if ( s_eqi ( command,'g') ) then call graph ( arrow,cp,delx,dely,dev,dudxn,dudyn,dvdxn,dvdyn,echo,eflag, & eflagu,etaref,filgrf,filtyp,icmax,icmin,icolor,icomp,iplot,isotri, & itable,iwork1,iwork2,jbound,jcmax,jcmin,lbar,line,maxbou,maxnp,maxnpe, & maxobj,maxsen,nbound,ncon,nelem,nflag,nflag0,node,np, & npe,nxskip,ny,nyskip,object,ovrlay,p,rho,rmach,s,s2,scalee, & scalen,scalev,show,smax,smin,srange,t,title,title2,u,v,x1max,x1min, & x2max,x2min,xc,xprof,xsiref,xsmax,xsmin,y1max,y1min,y2max,y2min, & yc,ysmax,ysmin) ! ! GRACE = : set the grace margin. ! else if ( s_eqi ( command(1:5),'grace') ) then if ( s_eqi ( command(1:6),'grace=') ) then read(command(7:),*,err = 90,end=90)grace else write(*,*)'Enter the grace margin:' read(*,*)grace write(17,*)grace if ( echo ) then write(*,*)grace end if end if write(*,*)' ' write(*,*)'DISPLAY - Note:' write(*,*)' The grace margin was set to GRACE = ',grace call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! HELLO: print program version, data, and maxima: ! else if ( s_eqi ( command,'hello') ) then call hello(maxbou,maxelm,maxnp,maxnpe,maxnx,maxny,maxobj,maxpar,maxsen) ! ! HELP: help ! else if ( s_eqi ( command,'h').or.s_eqi ( command,'help') ) then call help(echo) ! ! ICMAX = : set the maximum available color index. ! else if ( s_eqi ( command(1:6),'icmax=') ) then read(command(7:),*,err = 90,end=90)icmax if ( icmax > 255 ) then write(*,*)'ICMAX must be no more than 255' icmax = 255 end if write(*,*)'Maximum available color set to ',icmax ! ! ICMIN = : set the minimum available color index. ! else if ( s_eqi ( command(1:6),'icmin=') ) then read(command(7:),*,err = 90,end=90)icmin if ( icmin < 2 ) then write(*,*)'ICMIN must be at least 2.' icmin = 2 end if write(*,*)'Minimum available color set to ',icmin ! ! ICOMP = : set the sensitivity component. ! else if ( s_eqi ( command(1:5),'icomp') ) then if ( s_eqi ( command(1:6),'icomp=') ) then read(command(7:),*,err = 90,end=90)itemp else write(*,*)'Enter the sensitivity component:' read(*,*)itemp write(17,*)itemp if ( echo ) then write(*,*)itemp end if end if if ( 0 == itemp ) then write(*,*)'Primitive variables will be displayed.' else if ( 1 <= itemp.and.itemp <= npar ) then icomp = itemp write(*,*)'Sensitivity component set to ',icomp else if ( npar+1 <= itemp.and.itemp <= 2*npar ) then icomp = itemp write(*,*)'Finite difference component set to ',icomp-npar else write(*,*)'Display - Warning!' write(*,*)' Your command ICOMP = ',itemp write(*,*)' was ignored, because ICOMP must be' write(*,*)' 0 for the primitive variables,' write(*,*)' between 1 and ',npar,' for disc. sens.' write(*,*)' between ',npar+1,' and ',2*npar,' fin. dif.' end if ! ! INIT: set variables to initial (zero!) values. ! else if ( s_eqi ( command(1:4),'init') ) then call init ( arrow,cp,delx,dely,dev,dudxn,dudyn,dvdxn,dvdyn,echo,eflag, & eflagu,eqn,etaref,fildat,filgrf,filinp,filtyp,grace,icmax,icmin, & icolor,icomp,idata,ifile,iplot,iset,itable,iwrite,jcmax, & jcmin,labelx,labely,lbar,line,lppro,lptpro,lupro, & lutpro,lvpro,lvtpro,maxelm,maxnp,maxnpe,maxny,maxobj, & maxpar,maxsen,nbound,ncon,nelem,nflag,nflag0,node, & np,npe,nprof,nsen,nset,nxskip,ny,nyskip,object,ovrlay,p,rho, & rmach,scalee, & scalen,scalev,show,smax,smin,title,title2,u,v,x1max,x1min, & x2max,x2min,x4max,x4min,xc,xprof,xsiref,xsmax,xsmin,y1max, & y1min,y2max,y2min,y4max,y4min,yc,ysmax,ysmin) ! ! IWRITE = : set the debugging output level. ! else if ( s_eqi ( command(1:7),'iwrite=') ) then read(command(8:),*,err = 90,end=90)iwrite write(*,*)'Debugging level set to ',iwrite ! ! JCMAX = : set the maximum used color index. ! else if ( s_eqi ( command(1:6),'jcmax=') ) then read(command(7:),*,err = 90,end=90)jcmax if ( jcmax > 255 ) then jcmax = 255 write(*,*)'JCMAX must be no more than 255.' end if write(*,*)'Maximum used color set to ',jcmax ! ! JCMIN = : set the minimum used color index. ! else if ( s_eqi ( command(1:6),'jcmin=') ) then read(command(7:),*,err = 90,end=90)jcmin if ( jcmin < 2 ) then jcmin = 2 write(*,*)'JCMIN must be no less than 2.' end if write(*,*)'Minimum used color set to ',jcmin ! ! KV: show kinematic velocity vectors. ! else if ( s_eqi ( command,'kv') ) then show(8) = .not.show(8) if ( show(8) ) then write(*,*)'Kinematic velocities will be shown.' else write(*,*)'Kinematic velocities will NOT be shown.' end if ! ! KVMAG: show kinematic velocity magnitude contours. ! else if ( s_eqi ( command,'kvmag') ) then show(10) = .not.show(10) if ( show(10) ) then write(*,*)'Kinematic velocity magnitudes will be shown.' else write(*,*)'Kinematic velocity magnitudes will NOT be shown.' end if ! ! KVMAGC: show velocity magnitude color plots. ! else if ( s_eqi ( command,'kvmagc') ) then show(14) = .not.show(14) if ( show(14) ) then write(*,*)'Kinematic velocity magnitude colors will be shown.' else write(*,*)'Kinematic velocity magnitude colors will NOT be shown.' end if ! ! KVX: show X kinematic velocity contours. ! else if ( s_eqi ( command,'kvx') ) then show(15) = .not.show(15) if ( show(15) ) then write(*,*)'X kinematic velocity contours will be shown.' else write(*,*)'X kinematic velocity contours will NOT be shown.' end if ! ! KVXC: show X kinematic velocity color contours. ! else if ( s_eqi ( command,'kvxc') ) then show(17) = .not.show(17) if ( show(17) ) then write(*,*)'X kinematic velocity color contours will be shown.' else write(*,*)'X kinematic velocity color contours will NOT be shown.' end if ! ! KVY: show Y kinematic velocity contours. ! else if ( s_eqi ( command,'kvy') ) then show(16) = .not.show(16) if ( show(16) ) then write(*,*)'Y kinematic velocity contours will be shown.' else write(*,*)'Y kinematic velocity contours will NOT be shown.' end if ! ! KVYC: show Y kinematic velocity color contours. ! else if ( s_eqi ( command,'kvyc') ) then show(18) = .not.show(18) if ( show(18) ) then write(*,*)'Y kinematic velocity color contours will be shown.' else write(*,*)'Y kinematic velocity color contours will NOT be shown.' end if ! ! L: list current values ! else if ( (s_eqi ( command,'l')).or. & ( s_eqi (command(1:4),'list')) ) then call list(delx,dely,dev,echo,fildat,grace,icmax,icmin,icolor, & icomp,idata,ifile,iplot,iset,itable,iwrite,jbound, & maxbou,maxnp,maxobj,maxpar,maxsen,nbound,ncon,nelem,np,npar, & npe,nset,nx,nxskip,ny,nyskip,object,p,scalev,show,& title,title2,u,v,x2max,x2min,xmax,xmin,xprof,xsmax,xsmin,ymax, & ymin,y2max,y2min,ysmax,ysmin) ! ! LH: left half. ! else if ( s_eqi ( command,'lh') ) then xsmax = xsmin+0.5*(xsmax-xsmin) temp = 0.25*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! LINE: choose line type. ! else if ( s_eqi ( command,'line') ) then write(*,*)' ' write(*,*)' Number Linetype Name' write(*,*)' ' do i = 1,maxobj write(*,'(1x,i2,2x,i3,2x,a)')i,line(i),object(i) end do write(*,*)' ' write(*,*)'Enter an object number, and a line type.' write(*,*)'0 = Solid black, 1=dashed black,' write(*,*)'2 = Solid color, 3=dashed color.' read(*,*,end = 90,err=90)itemp,jtemp write(17,*)itemp,jtemp if ( echo ) then write(*,*)itemp,jtemp end if if ( 1 <= itemp.and.itemp <= maxobj ) then line(itemp) = jtemp else write(*,*)'Your object number was out of bounds.' end if ! ! MACH: draw mach contour lines. ! else if ( s_eqi ( command,'mach') ) then show(30) = .not.show(30) if ( show(30) ) then write(*,*)'Mach contours will be plotted.' else write(*,*)'Mach contours will NOT be plotted.' end if ! ! MACHC: draw density colors. ! else if ( s_eqi ( command,'machc') ) then show(33) = .not.show(33) if ( show(33) ) then write(*,*)'Mach colors will be plotted.' else write(*,*)'Mach colors will NOT be plotted.' end if ! ! MC: middle center quarter. ! else if ( s_eqi ( command,'mc') ) then temp = 0.25*(xsmax-xsmin) xsmin = xsmin+temp xsmax = xsmax-temp temp = 0.25*(ysmax-ysmin) ysmax = ysmax-temp ysmin = ysmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! MH: middle half. ! else if ( s_eqi ( command,'mh') ) then temp = 0.25*(ysmax-ysmin) ysmax = ysmax-temp ysmin = ysmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! ML: middle left quarter. ! else if ( s_eqi ( command,'ml') ) then xsmax = xsmin+0.5*(xsmax-xsmin) temp = 0.25*(ysmax-ysmin) ysmax = ysmax-temp ysmin = ysmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! MR: middle right quarter. ! else if ( s_eqi ( command,'mr') ) then xsmin = xsmin+0.5*(xsmax-xsmin) temp = 0.25*(ysmax-ysmin) ysmax = ysmax-temp ysmin = ysmin+temp call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! MV: show kinematic velocity vectors. ! else if ( s_eqi ( command,'mv') ) then show(35) = .not.show(35) if ( show(35) ) then write(*,*)'Mass velocities will be shown.' else write(*,*)'Mass velocities will NOT be shown.' end if ! ! MVMAG: show mass velocity magnitude contours. ! else if ( s_eqi ( command,'mvmag') ) then show(36) = .not.show(36) if ( show(36) ) then write(*,*)'Mass velocity magnitudes will be shown.' else write(*,*)'Mass velocity magnitudes will NOT be shown.' end if ! ! MVMAGC: show mass magnitude color plots. ! else if ( s_eqi ( command,'mvmagc') ) then show(37) = .not.show(37) if ( show(37) ) then write(*,*)'Mass velocity magnitude colors will be shown.' else write(*,*)'Mass velocity magnitude colors will NOT be shown.' end if ! ! N: show nodes ! else if ( s_eqi ( command,'n') ) then show(4) = .not.show(4) if ( show(4) ) then write(*,*)'Nodes will be shown.' else write(*,*)'Nodes will NOT be shown.' end if ! ! NCON = : Set the number of contour lines. ! else if ( s_eqi ( command(1:4),'ncon') ) then if ( s_eqi ( command(1:5),'ncon=') ) then read(command(6:),*,err = 90,end=90)ncon write(*,*)'Number of contour lines set to ',ncon else write(*,*)'Enter number of contour lines.' read(*,*)ncon write(17,*)ncon if ( echo ) then write(*,*)ncon end if end if ! ! NN: Show node numbers ! else if ( s_eqi ( command,'nn') ) then show(27) = .not.show(27) if ( show(27) ) then write(*,*)'Node numbers will be shown.' else write(*,*)'Node numbers will NOT be shown.' end if ! ! NOFRAME: no frame, please ! else if ( s_eqi ( command,'noframe') ) then show(3) = .false. write(*,*)'No frame will be shown around the picture.' ! ! NXSKIP = : skip value for column nodes. ! else if ( s_eqi ( command(1:7),'nxskip=') ) then read(command(8:),*,err = 90,end=90)itemp if ( itemp <= 0 ) then write(*,*)'Error! NXSKIP must be 1 or more!' else nxskip = itemp write(*,*)'Skip value for column nodes is ',nxskip end if ! ! NYSKIP = : skip value for row nodes. ! else if ( s_eqi ( command(1:7),'nyskip=') ) then read(command(8:),*,err = 90,end=90)itemp if ( itemp <= 0 ) then write(*,*)'Error! NYSKIP must be 1 or more!' else nyskip = itemp write(*,*)'Skip value for row nodes is ',nyskip end if ! ! OVERLAY: Switch the overlay value. ! else if ( s_eqi ( command,'overlay') ) then ovrlay = .not.ovrlay if ( ovrlay ) then write(*,*)'Plots will be overlayed until next OVERLAY.' else write(*,*)'This overlay plot is done.' call newfrm end if ! ! P: show pressure contours. ! else if ( s_eqi ( command,'p') ) then show(5) = .not.show(5) if ( show(5) ) then write(*,*)'Pressures will be shown.' else write(*,*)'Pressures will NOT be shown.' end if ! ! PC: show pressure color plots. ! else if ( s_eqi ( command,'pc') ) then show(12) = .not.show(12) if ( show(12) ) then write(*,*)'Pressure colors will be shown.' else write(*,*)'Pressure colors will NOT be shown.' end if ! ! PPRO: draw P profile ! else if ( s_eqi ( command,'ppro') ) then lppro = .not.lppro if ( lppro ) then write(*,*)'Pressure will be shown in profile.' else write(*,*)'Pressure will NOT be shown in profile.' end if ! ! PROGRAF: show a graph of the profile data. ! else if ( s_eqi ( command,'prograf') ) then call pgraf(dev,echo,filgrf,icmax,icmin,icolor,iplot,itable, & labelx,labely,lppro,lptpro,lupro,lutpro,lvpro,lvtpro,maxnp,maxobj, & maxsen,np,nprof,ny,ovrlay,p,ptar,show,title,title2,u,utar,v, & vtar,x1max,x1min,x2max,x2min,x4max,x4min,y1max,y1min,y2max, & y2min,y4max,y4min,yc) ! ! PROF: show profile line. ! else if ( s_eqi ( command(1:4),'prof') ) then show(19) = .not.show(19) if ( show(19) ) then write(*,*)'The profile line will be shown.' else write(*,*)'The profile line will NOT be shown.' end if ! ! PTPRO: draw P target profile ! else if ( s_eqi ( command,'ptpro') ) then lptpro = .not.lptpro if ( lptpro ) then write(*,*)'Target pressure will be shown in profile.' else write(*,*)'Target pressure will NOT be shown in profile.' end if ! ! Q: quit. ! QU: QUIT NOW! ! else if ( s_eqi ( command(1:1), 'Q' ) ) then if ( s_eqi ( command(2:2), 'U' ) ) then command = 'Y' else write ( *, * ) 'Enter "Y" to confirm you want to quit!' read ( *, '(a)' ) command write ( 17, '(a)' ) trim ( command ) if ( echo ) then write(*,'(a)') trim ( command ) end if end if if ( s_eqi ( command, 'Y' ) ) then call quit ( dev, ifile, iplot ) end if ! ! RH: right half. ! else if ( s_eqi ( command,'rh') ) then temp = 0.50*(xsmax-xsmin) xsmin = xsmin+temp call pltbox ( grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! RHO: draw density contour lines. ! else if ( s_eqi ( command,'rho') ) then show(31) = .not.show(31) if ( show(31) ) then write(*,*)'Rho contours will be plotted.' else write(*,*)'Rho contours will NOT be plotted.' end if ! ! RHOC: draw density colors. ! else if ( s_eqi ( command,'rhoc') ) then show(34) = .not.show(34) if ( show(34) ) then write(*,*)'Rho colors will be plotted.' else write(*,*)'Rho colors will NOT be plotted.' end if ! ! S: show stream lines ! else if ( s_eqi ( command,'s').or. s_eqi ( command(1:6),'stream') ) then show(6) = .not.show(6) if ( show(6) ) then write(*,*)'Stream lines will be shown.' else write(*,*)'Stream lines will NOT be shown.' end if ! ! SC: show stream colors ! else if ( s_eqi ( command,'sc') ) then show(22) = .not.show(22) if ( show(22) ) then write(*,*)'Stream colors will be shown.' else write(*,*)'Stream colors will NOT be shown.' end if ! ! SCALEE = : set the element number scale factor. ! else if ( s_eqi ( command(1:7),'scalee=') ) then read ( command(8:),*,err = 90,end=90)scalee write(*,*)'Element number scale factor set to ',scalee ! ! SCALEN = : set the node number scale factor. ! else if ( s_eqi ( command(1:7),'scalen=') ) then read(command(8:),*,err = 90,end=90)scalen write(*,*)'Node number scale factor set to ',scalen ! ! SCALEV = : set the velocity vector scale. ! else if ( s_eqi ( command(1:7),'scalev=') ) then read(command(8:),*,err = 90,end=90)scalev write(*,*)'Velocity vector scale factor set to ',scalev ! ! SMAX = : enter contour maximum value. ! else if ( s_eqi ( command(1:5),'smax=') ) then read(command(6:),*,err = 90,end=90)smax else if ( s_eqi ( command(1:4),'smax') ) then write(*,*)'Enter the contour maximum range:' read(*,*,end = 80)smax write(17,*)smax if ( echo ) then write(*,*)smax end if ! ! SMIN = : enter contour minimum value. ! else if ( s_eqi ( command(1:5),'smin=') ) then read(command(6:),*,err = 90,end=90)smin else if ( s_eqi ( command(1:4),'smin') ) then write(*,*)'Enter the contour minimum range:' read(*,*,end = 80)smin write(17,*)smin if ( echo ) then write(*,*)smin end if ! ! TC: top center quarter. ! else if ( s_eqi ( command,'tc') ) then temp = 0.25*(xsmax-xsmin) xsmin = xsmin+temp xsmax = xsmax-temp ysmin = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! TH: top half ! else if ( s_eqi ( command,'th') ) then ysmin = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! TITLE2 = : enter subtitle ! else if ( s_eqi ( command(1:7),'title2=') ) then title2 = command(8:) else if ( s_eqi ( command(1:6),'title2') ) then write(*,*)'Enter the subtitle:' read(*,'(a)',end = 80)title2 write(17,'(a)') trim ( title2 ) if ( echo ) then write(*,'(a)') trim ( title2 ) end if ! ! TITLE = : enter title ! else if ( s_eqi ( command(1:6),'title=') ) then title = command(7:) else if ( s_eqi ( command(1:5),'title') ) then write(*,*)'Enter the plot title:' read(*,'(a)',end = 80)title write(17,'(a)') trim ( title ) if ( echo ) then write(*,'(a)') trim ( title ) end if ! ! TL: top left quarter. ! else if ( s_eqi ( command,'tl') ) then xsmax = xsmin+0.5*(xsmax-xsmin) ysmin = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! TR: top right quarter. ! else if ( s_eqi ( command,'tr') ) then xsmin = xsmin+0.5*(xsmax-xsmin) ysmin = ysmin+0.5*(ysmax-ysmin) call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax,xsmin, & y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! UPRO: draw U profile ! else if ( s_eqi ( command,'upro') ) then lupro = .not.lupro if ( lupro ) then write(*,*)'Horizontal velocity will be shown in profile.' else write(*,*)'Horizontal velocity will NOT be shown in profile.' end if ! ! UTPRO: draw U target profile ! else if ( s_eqi ( command,'utpro') ) then lutpro = .not.lutpro if ( lutpro ) then write(*,*)'Target horizontal velocity will be shown in profile.' else write(*,*)'Target horizontal velocity will NOT be shown in profile.' end if ! ! UV: show unit velocity direction vectors. ! else if ( s_eqi ( command,'uv') ) then show(9) = .not.show(9) if ( show(9) ) then write(*,*)'Unit velocity vectors will be shown.' else write(*,*)'Unit velocity vectors will NOT be shown.' end if ! ! VE: set visible elements. ! else if ( s_eqi ( command,'ve') ) then call element_visibility(eflagu,nelem) write(*,*)'Do you want to adjust the data limits' write(*,*)'to focus on the visible elements?' read(*,'(a1)')isay write(17,'(a1)')isay if ( echo ) then write(*,'(a1)')isay end if if ( s_eqi ( isay,'y') ) then none = .true. do i = 1,nelem if ( eflagu(i) ) then do j = 1,npe k = node(j,i) if ( none ) then xsmin = xc(k) xsmax = xc(k) ysmin = yc(k) ysmax = yc(k) none = .false. else xsmin = min(xsmin,xc(k)) xsmax = max(xsmax,xc(k)) ysmin = min(ysmin,yc(k)) ysmax = max(ysmax,yc(k)) end if end do end if end do if ( none ) then write(*,*)' ' write(*,*)'DISPLAY - Warning!' write(*,*)' No visible elements were found!' else call pltbox(grace,srange,x1max,x1min,x2max,x2min,xsmax, & xsmin,y1max,y1min,y2max,y2min,ysmax,ysmin) end if end if ! ! VN: set visible nodes. ! else if ( s_eqi ( command,'vn') ) then call viznod ( echo, np, nflag0 ) ! ! VND: set visible nodes by distance. ! else if ( s_eqi ( command,'vnd') ) then call viznd ( echo, np, nflag0, xc, yc ) ! ! VORT: show vorticity ! else if ( s_eqi ( command,'vort') ) then show(11) = .not.show(11) if ( show(11) ) then write(*,*)'Vorticity contours will be shown.' else write(*,*)'Vorticity contours will NOT be shown.' end if ! ! VORTC: show vorticity color plots. ! else if ( s_eqi ( command,'vortc') ) then show(13) = .not.show(13) if ( show(13) ) then write(*,*)'Vorticity colors will be shown.' else write(*,*)'Vorticity colors will NOT be shown.' end if ! ! VPRO: draw V profile ! else if ( s_eqi ( command,'vpro') ) then lvpro = .not.lvpro if ( lvpro ) then write(*,*)'Vertical velocity will be shown in profile.' else write(*,*)'Vertical velocity will NOT be shown in profile.' end if ! ! VTPRO: draw V target profile ! else if ( s_eqi ( command,'vtpro') ) then lvtpro = .not.lvtpro if ( lvtpro ) then write(*,*)'Target vertical velocity will be shown in profile.' else write(*,*)'Target vertical velocity will NOT be shown in profile.' end if ! ! X: set the data window. ! else if ( s_eqi ( command,'x') ) then call get_win(echo,grace,srange,xmax,xmin,x1max,x1min,x2max, & x2min,xsmax,xsmin,ymax,ymin,y1max,y1min,y2max,y2min,ysmax,ysmin) ! ! XC: show X coordinate contours. ! else if ( s_eqi ( command,'xc') ) then show(25) = .not.show(25) if ( show(25) ) then write(*,*)'X coordinate contours will be shown.' else write(*,*)'X coordinate contours will NOT be shown.' end if ! ! XCC: show X coordinate color plots. ! else if ( s_eqi ( command,'xcc') ) then show(23) = .not.show(23) if ( show(23) ) then write(*,*)'X coordinate colors will be shown.' else write(*,*)'X coordinate colors will NOT be shown.' end if ! ! XPROF = : reset location of profile line (til next step is read!). ! else if ( s_eqi ( command,'xprof=') ) then read(command(7:),*,err = 90,end=90)xprof write(*,*)'X coordinate of profile reset to ',xprof ! ! YC: show Y coordinate contours. ! else if ( s_eqi ( command,'yc') ) then show(26) = .not.show(26) if ( show(26) ) then write(*,*)'Y coordinate contours will be shown.' else write(*,*)'Y coordinate contours will NOT be shown.' end if ! ! YCC: show Y coordinate color plots. ! else if ( s_eqi ( command,'ycc') ) then show(24) = .not.show(24) if ( show(24) ) then write(*,*)'Y coordinate colors will be shown.' else write(*,*)'Y coordinate colors will NOT be shown.' end if ! ! #: a comment ! else if ( command(1:1) == '#' ) then ! ! Unrecognized command. ! else write(*,*)' ' write(*,*)'DISPLAY - Warning!' write(*,*)' DISPLAY did not recognize your command:' write(*,'(2x,''"'',a,''"'')') trim ( command ) end if ! ! Go back, and get next command. ! go to 10 80 continue write(*,*)' ' write(*,*)'DISPLAY - Fatal error!' write(*,*)' Error or end of file reading user input!' call quit(dev,ifile,iplot) 90 continue write(*,*)' ' write(*,*)'DISPLAY - Serious error!' write(*,*)' Error reading user input!' command = ' ' go to 10 100 continue write(*,*)' ' write(*,*)'DISPLAY - Serious error!' write(*,*)' The input file you named could not be opened!' command = ' ' go to 10 110 continue write(*,*)' ' write(*,*)'DISPLAY - Serious error!' write(*,*)' End of file, or error occurred, while trying' write(*,*)' to read the input file you named.' command = ' ' go to 10 end subroutine adjeff ( delx,dely,eflag,eqn,fildat,grace,ifile,isotri,jbound, & jfile,maxbou,maxelm,maxnp,maxnpe,nbound,nelem,nflag,nflag0,node,np,npe, & nset,nx,ny,p,rho,srange,u,v,x1max,x1min,x2max,x2min,xc,xmax,xmin,xsmax, & xsmin,xtmax,xtmin,y1max,y1min,y2max,y2min,yc,ymax,ymin,ysmax,ysmin, & ytmax,ytmin) ! !*********************************************************************** ! !! ADJEFF reads the node data for a Jeff plot. ! ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real DELX, the X spacing between nodes. In some cases, ! this spacing is modified to create isoparametric elements. ! ! Input, real DELY, the Y spacing between nodes. In some cases, ! this spacing is modified to create isoparametric elements. ! ! Output, logical EFLAG(MAXELM). ! EFLAG is used to "flag" which elements are to be displayed. ! If EFLAG(I) is TRUE, then element I should be displayed. ! ! Output, character ( len = 2 ) EQN(3,NP). ! EQN records the "type" of each equation that will be generated, and ! which is associated with an unknown. Note that most boundary ! conditions do not result in an equation. The current values are: ! ! 'U' The horizontal momentum equation. ! 'UW' The condition U = 0 applied at a node on a fixed wall. ! 'V' The vertical momentum equation. ! 'VW' The condition V = 0 applied at a node on a fixed wall. ! 'P' The continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! ! Input, character ( len = 80 ) FILDAT. ! The name of the data file to be read in, which contains ! the information defining the mesh and the physical ! parameters. ! ! Input, real GRACE. ! The size of the "grace" margin on the plot. ! ! IFILE Input/output, integer IFILE. ! Records the status of the data file whose name is FILDAT. ! ! ISOTRI Output, integer ISOTRI(MAXELM). ! ! 0, if element I is not isoparametric. ! 1 or 2, if element I is isoparametric. ! ! JBOUND Output, integer JBOUND(5,MAXBOU) ! ! For each line segment of the boundary: ! ! JBOUND(1,I) contains the element number; ! ! JBOUND(2,I) contains the local node number of one corner ! of the element, which forms the edge; ! ! JBOUND(2,I) contains the "next" node along the edge. ! If the element is linear, this is the other corner node. ! If the element is quadratic, this is the midside node along ! the edge. ! ! JBOUND(4,I) contains the "next" node along the edge. ! If the element is linear, this is 0. ! If the element is quadratic, this is the other corner node ! along the edge. ! ! JBOUND(5,I) contains: ! 0 if the boundary is a wall (U = V=0); ! 1 if the boundary is open. ! ! JFILE Output, integer JFILE, an error indicator. ! JFILE is 2 if the file was successfully opened and read. ! (If not, then RdNod actually halts execution!) ! ! MAXBOU Input, integer MAXBOU. ! The amount of storage available for the IBOUND array. ! ! MAXELM Input, integer MAXELM. ! The maximum number of elements which the program can ! handle. ! ! MAXNP Input, integer MAXNP, the maximum number of nodes. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! MAXNY Input, integer MAXNY, the maximum allowed value of NY. ! ! NBOUND Output, integer NBOUND. ! The number of points (XBOUND(I),YBOUND(I)) used to ! define the boundary. ! ! NELEM Output, integer NELEM. ! The number of elements. ! ! NFLAG Output, logical NFLAG(MAXNP). ! ! NFLAG is used to "flag" which nodes are active, ! that is, to be displayed, and which not, in the graph. ! ! NODE Output, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Output, integer NP, the number of nodes. ! ! NPE Output, integer NPE. ! NPE is the number of nodes per element, which should be ! 3 for linear elements and 6 for quadratics. ! ! NSET Output, integer NSET. ! The number of sets of data contained in the data file. ! ! NX Output, integer NX. ! Determines the number of nodes and elements in the X ! direction. There will be 2*NX-1 nodes, 2*NX-2 elements. ! ! NY Output, integer NY. ! Determines the number of nodes and elements in the Y ! direction. There will be 2*NY-1 nodes, 2*NY-2 elements. ! ! P Output, real P(MAXNP,0:MAXPAR). ! ! P(I,0) is the pressure at node I. ! ! P(I,J) is the sensitivity of the pressure with respect ! to parameter J. ! ! SRANGE Output, real SRANGE. ! The maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! U Output, real U(MAXNP,MAXPAR). ! ! U(I,0) is the horizontal fluid velocity at node I. ! ! U(I,J) is the sensitivity of the horizontal velocity with ! respect to parameter J. ! ! V Output, real V(MAXNP,MAXPAR). ! ! V(I,0) is the vertical fluid velocity at node I. ! ! V(I,J) is the sensitivity of the vertical velocity with ! respect to parameter J. ! ! X1MAX, ! X1MIN Output, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! X2MAX, ! X2MIN Output, real X2MAX, X2MIN, the maximum and minimum X ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! XC Output, real XC(MAXNP). ! XC contains the X coordinates of the nodes. ! ! XMAX Output, real XMAX. ! The maximum X coordinate of all the nodes. ! The maximum entry in the XC array. ! ! XMIN Output, real XMIN. ! The minimum X coordinate of all the nodes. ! The minimum entry in the XC array. ! ! XSMAX Output, real XSMAX. ! The maximum X coordinate of the data to be displayed. ! XSMAX defaults to XMAX, but can be made smaller to ! focus on a portion of the region. ! ! XSMIN Output, real XSMIN. ! The minimum X coordinate of the data to be displayed. ! XSMIN defaults to XMIN, but can be made larger to ! focus on a portion of the region. ! ! Y1MAX, ! Y1MIN Output, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! ! Y2MAX, ! Y2MIN Output, real Y2MAX, Y2MIN, the maximum and minimum Y ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! YC Output, real YC(MAXNP). ! The Y coordinates of the nodes. ! ! YMAX Output, real YMAX. ! The maximum Y coordinate of all the nodes. ! The maximum value attained by the YC array. ! ! YMIN Output, real YMIN. ! The minimum Y coordinate of all the nodes. ! The minimum value attained by the YC array. ! ! YSMAX Output, real YSMAX. ! The maximum Y coordinate of the data to be displayed. ! YSMAX defaults to YMAX, but can be made smaller to ! focus on a portion of the region. ! ! YSMIN Output, real YSMIN. ! The minimum Y coordinate of the data to be displayed. ! YSMIN defaults to YMIN, but can be made larger to ! focus on a portion of the region. ! integer maxbou integer maxelm integer maxnp integer maxnpe ! real delx real dely logical eflag(maxelm) character ( len = 2 ) eqn(3,maxnp) character ( len = 80 ) fildat real grace integer i integer ifile logical inside integer isotri(maxelm) integer jbound(5,maxbou) integer jfile integer nbound integer nelem logical nflag(maxnp) logical nflag0(maxnp) integer node(maxnpe,maxelm) integer np integer npe integer nset integer nx integer ny real p(maxnp) real rho(maxnp) real srange real u(maxnp) real v(maxnp) real xc(maxnp) real xmax real xmin real x1max real x1min real x2max real x2min real xsmax real xsmin real xtmax real xtmin real yc(maxnp) real ymax real ymin real y1max real y1min real y2max real y2min real ysmax real ysmin real ytmax real ytmin ! ! Read nodal data. ! call rdnod ( fildat, jfile, np, p, rho, u, v, xc, yc ) nset = 1 ifile = 2 ! ! The values for NX, and NY are just dummy values for now. ! nx = 1 ny = 1 ! ! We assume all triangles are nonisoperimetric. ! isotri(1:nelem) = 0 ! ! Initially, all elements will be visible. ! do i = 1,nelem eflag(i) = .true. end do ! do i = 1,np nflag(i) = .true. nflag0(i) = .true. end do ! ! We give default values to the equations. ! do i = 1,np eqn(1,i) = 'U' eqn(2,i) = 'V' eqn(3,i) = 'P' end do write(*,*)' ' write(*,*)'AdJeff - Note:' write(*,*)' The node data has been read.' ! ! Check element orientation. ! call element_check(maxnpe,nelem,node,np,npe,xc,yc) ! ! Determine the region size. ! call region_size ( delx,dely,grace,nelem,nflag,np,srange, & x1max,x1min,x2max,x2min,xc,xmax,xmin,xsmax,xsmin,xtmax,xtmin,y1max, & y1min,y2max,y2min,yc,ymax,ymin,ysmax,ysmin,ytmax,ytmin) ! ! Set the location of boundary edges. ! call boundary_set(eqn,jbound,maxbou,maxnpe,nbound,nelem,node,np,npe) return end subroutine advanc ( delx,dely,echo,eflag,eqn,grace,ifile,iset,isotri, & jbound,maxbou,maxelm,maxnp,maxnpe,maxny,maxpar,maxsen,nbound, & nelem,nflag,nflag0,node,np,npar,npe,nprof,nsen,nset,nx,ny,p, & para,ptar,rho,srange,u,utar,v,vtar,x1max,x1min,x2max,x2min,xc, & xmax,xmin,xsmax,xsmin,xtmax,xtmin,y1max,y1min,y2max,y2min,yc, & ymax,ymin,ysmax,ysmin,ytmax,ytmin) ! !*********************************************************************** ! !! ADVANC responds to a user's "A" command by asking which step we are ! to advance to, and then reading the information for that step. ! ! ! Modified: ! ! 07 January 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real DELX, the X spacing between nodes. In some cases, ! this spacing is modified to create isoparametric elements. ! ! Output, real DELY, the Y spacing between nodes. In some cases, ! this spacing is modified to create isoparametric elements. ! ! Output, logical EFLAG(MAXELM). ! EFLAG is used to "flag" which elements are to be displayed. ! If EFLAG(I) is TRUE, then element I should be displayed. ! ! Output, character ( len = 2 ) EQN(3,NP). ! EQN records the "type" of each equation that will be generated, and ! which is associated with an unknown. Note that most boundary ! conditions do not result in an equation. The current values are: ! ! 'U' The horizontal momentum equation. ! 'UW' The condition U = 0 applied at a node on a fixed wall. ! 'V' The vertical momentum equation. ! 'VW' The condition V = 0 applied at a node on a fixed wall. ! 'P' The continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! ! Input, real GRACE. ! The size of the "grace" margin on the plot. ! ! IFILE Input/output, integer IFILE. ! Records the status of the data file whose name is FILDAT. ! ! -2, an error occurred while reading from the file. ! -1, the file could not be opened. ! 0, no file is currently open. ! 1, a file has been opened, but not read from. ! 2, data has been read from a file. ! ! ISET Output, integer ISET. ! The data set being examined from the file. If no file ! is open, or if no data set has been read, then ISET is ! zero. ! ! ISOTRI Output, integer ISOTRI(MAXELM). ! ! 0, if element I is not isoparametric. ! 1 or 2, if element I is isoparametric. ! ! JBOUND Output, integer JBOUND(5,MAXBOU) ! ! For each line segment of the boundary: ! ! JBOUND(1,I) contains the element number; ! ! JBOUND(2,I) contains the local node number of one corner ! of the element, which forms the edge; ! ! JBOUND(2,I) contains the "next" node along the edge. ! If the element is linear, this is the other corner node. ! If the element is quadratic, this is the midside node along ! the edge. ! ! JBOUND(4,I) contains the "next" node along the edge. ! If the element is linear, this is 0. ! If the element is quadratic, this is the other corner node ! along the edge. ! ! JBOUND(5,I) contains: ! 0 if the boundary is a wall (U = V=0); ! 1 if the boundary is open. ! ! MAXBOU Input, integer MAXBOU. ! The amount of storage available for the IBOUND array. ! ! MAXELM Input, integer MAXELM. ! The maximum number of elements which the program can ! handle. ! ! MAXNP Input, integer MAXNP, the maximum number of nodes. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! MAXPAR Input, integer MAXPAR. ! The maximum number of parameters the program can handle. ! ! NBOUND Output, integer NBOUND. ! The number of points (XBOUND(I),YBOUND(I)) used to ! define the boundary. ! ! NELEM Output, integer NELEM. ! The number of elements. ! ! NFLAG Output, logical NFLAG(MAXNP). ! ! NFLAG is used to "flag" which nodes are active, ! that is, to be displayed, and which not, in the graph. ! ! NODE Output, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Output, integer NP, the number of nodes. ! ! NPAR Output, integer NPAR. ! The number of parameters. ! ! NPE Output, integer NPE. ! NPE is the number of nodes per element, which should be ! 3 for linear elements and 6 for quadratics. ! ! NPROF Output, integer NPROF(MY). ! Records the indices of the nodes that lie along the profile ! line. ! ! NSET Output, integer NSET. ! The number of sets of data contained in the data file. ! ! NX Output, integer NX. ! Determines the number of nodes and elements in the X ! direction. There will be 2*NX-1 nodes, 2*NX-2 elements. ! ! NY Output, integer NY. ! Determines the number of nodes and elements in the Y ! direction. There will be 2*NY-1 nodes, 2*NY-2 elements. ! ! P Output, real P(MAXNP,0:MAXPAR). ! ! P(I,0) is the pressure at node I. ! ! P(I,J) is the sensitivity of the pressure with respect ! to parameter J. ! ! PARA Output, real PARA(MAXPAR). ! The value of the parameters. ! ! PTAR Output, real PTAR(MAXNP) ! The pressure field associated with the target solution, at ! node I. ! ! SRANGE Output, real SRANGE. ! The maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! U Output, real U(MAXNP,MAXPAR). ! ! U(I,0) is the horizontal fluid velocity at node I. ! ! U(I,J) is the sensitivity of the horizontal velocity with ! respect to parameter J. ! ! UTAR Output, real UTAR(MAXNP) ! The horizontal velocity field associated with the target ! solution, at node I. ! ! V Output, real V(MAXNP,MAXPAR). ! ! V(I,0) is the vertical fluid velocity at node I. ! ! V(I,J) is the sensitivity of the vertical velocity with ! respect to parameter J. ! ! VTAR Output, real VTAR(MAXNP) ! The vertical velocity field associated with the target solution, ! at node I. ! ! X1MAX, ! X1MIN Output, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! X2MAX, ! X2MIN Output, real X2MAX, X2MIN, the maximum and minimum X ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! XC Output, real XC(MAXNP). ! The X coordinates of the nodes. ! ! XMAX Output, real XMAX. ! The maximum X coordinate of all the nodes. ! The maximum entry in the XC array. ! ! XMIN Output, real XMIN. ! The minimum X coordinate of all the nodes. ! The minimum entry in the XC array. ! ! XSMAX Output, real XSMAX. ! The maximum X coordinate of the data to be displayed. ! XSMAX defaults to XMAX, but can be made smaller to ! focus on a portion of the region. ! ! XSMIN Output, real XSMIN. ! The minimum X coordinate of the data to be displayed. ! XSMIN defaults to XMIN, but can be made larger to ! focus on a portion of the region. ! ! Y1MAX, ! Y1MIN Output, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! ! Y2MAX, ! Y2MIN Output, real Y2MAX, Y2MIN, the maximum and minimum Y ! coordinates that should be used for plotting. No plotting ! commands should exceed these values. This is where the ! "frame" might be drawn. ! ! YC Output, real YC(MAXNP). ! The Y coordinates of the nodes. ! ! YMAX Output, real YMAX. ! The maximum Y coordinate of all the nodes. ! The maximum value attained by the YC array. ! ! YMIN Output, real YMIN. ! The minimum Y coordinate of all the nodes. ! The minimum value attained by the YC array. ! ! YSMAX Output, real YSMAX. ! The maximum Y coordinate of the data to be displayed. ! YSMAX defaults to YMAX, but can be made smaller to ! focus on a portion of the region. ! ! YSMIN Output, real YSMIN. ! The minimum Y coordinate of the data to be displayed. ! YSMIN defaults to YMIN, but can be made larger to ! focus on a portion of the region. ! integer maxbou integer maxelm integer maxnp integer maxnpe integer maxny integer maxpar integer maxsen ! character ( len = 6 ) chrint real delx real dely logical echo logical eflag(maxelm) character ( len = 2 ) eqn(3,maxnp) real grace integer i integer ierror integer ifile logical inside integer iset integer isotri(maxelm) integer iwant integer jbound(5,maxbou) integer nbound integer nelem logical nflag(maxnp) logical nflag0(maxnp) integer node(maxnpe,maxelm) integer np integer npar integer npe integer nprof(2*maxny-1) integer nsen integer nset integer nx integer ny real p(maxnp,0:maxsen) real para(maxpar) real ptar(maxnp) real rho(maxnp) real srange character ( len = 80 ) string real u(maxnp,0:maxsen) real utar(maxnp) real v(maxnp,0:maxsen) real vtar(maxnp) real x1max real x1min real x2max real x2min real xc(maxnp) real xmax real xmin real xprof real xsmax real xsmin real xtmax real xtmin real y1max real y1min real y2max real y2min real yc(maxnp) real ymax real ymin real ysmax real ysmin real ytmax real ytmin ! ! ADVANC expects to use 6 node quadrilaterals. ! npe = 6 ! if ( nset <= 0 ) then write(*,*)' ' write(*,*)'Advanc - Error!' write(*,*)' Please use the "DAT" command first, to define' write(*,*)' a data file to read!' return end if ! ! Get the desired step number. ! 10 continue write(*,*)' ' string = 'There are a total of '//chrint(nset)//' steps in the data file.' call s_blanks_delete ( string) write(*,'(1x,a)') write(*,*)'What step would you like to see?' read(*,*,err = 60,end=60)iwant write(17,*)iwant if ( echo ) then write(*,*)iwant end if if ( iwant <= 0 ) then write(*,*)' ' write(*,*)'Sorry, the step must be positive!' go to 10 else if ( iwant > nset ) then write(*,*)' ' string = 'Sorry, the step must be no greater than '// chrint(nset) call s_blanks_delete ( string) write(*,'(1x,a)')string go to 10 end if ! ! Advance to the desired step. ! ! If this REWIND doesn't work, you may have to close and then reopen. ! rewind 2 do iset = 0, iwant-1 call plot_file_read ( eqn,ierror,ifile,iset,isotri,maxelm,maxnp, & maxnpe,maxny,maxpar,maxsen,nelem,node,np,npar,npe,nprof, & nsen,nx,ny,p,para,ptar,rho,u,utar,v,vtar,xc,xprof,yc) if ( ierror /= 0 ) then write(*,*)' ' write(*,*)'ADVANC - Serious error!' write(*,*)' Error occurred while trying to reach a desired' write(*,*)' step. We were reading information for step' write(*,*)iset return end if end do nflag(1:np) = .true. nflag0(1:np) = .true. eflag(1:nelem) = .true. write(*,*)' ' write(*,*)'ADVANC - Note:' write(*,*)' Read data for step ',iwant ! ! Check the orientation of the elements. ! call element_check ( maxnpe, nelem, node, np, npe, xc, yc ) ! ! Determine the region size. ! call region_size ( delx,dely,grace,nelem,nflag,np,srange, & x1max,x1min,x2max,x2min,xc,xmax,xmin,xsmax,xsmin,xtmax,xtmin,y1max, & y1min,y2max,y2min,yc,ymax,ymin,ysmax,ysmin,ytmax,ytmin ) ! ! Compute the location of boundary edges. ! call boundary_set ( eqn, jbound, maxbou, maxnpe, nbound, nelem, node, & np, npe ) return 60 continue write(*,*)' ' write(*,*)'ADVANC - Serious error!' write(*,*)' Error while reading user input.' write(*,*)' ' return end subroutine arrlin ( xstart, ystart, xtip, ytip ) ! !*********************************************************************** ! !! ARRLIN can be used to make a line drawing of an arrow at any point ! on a graph. ! ! The arrow will stretch between two user specified points. ! ! The "head" of the arrow may be fatter or thinner than expected ! if the X and Y scales of the graph are not in the same ! proportions. ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! XSTART, ! YSTART Input, real XSTART, YSTART, the starting point for the ! arrow. ! ! XTIP, ! YTIP Input, real XTIP, YTIP, the end point for the arrow. ! ! ! left ! |\ ! | \ ! start ------- base tip ! | / ! |/ ! rite ! real pi parameter (pi = 3.1415926) ! real alpha real del real dist real theta real xbase real xleft real xrite real xstart real xtip real ybase real yleft real yrite real ystart real ytip ! ! if ( xstart == xtip.and.ystart == ytip)return theta = 0.5*pi-atan2(2.0,1.0) dist = sqrt((xtip-xstart)**2+(ytip-ystart)**2) alpha = atan2(ytip-ystart,xtip-xstart) del = sqrt(5.0)*dist/3.0 call movcgm(xstart,ystart) xbase = xstart+dist*cos(alpha)*2.0/3.0 ybase = ystart+dist*sin(alpha)*2.0/3.0 call drwcgm(xbase,ybase) xleft = xstart+del*cos(alpha-theta) yleft = ystart+del*sin(alpha-theta) call drwcgm(xleft,yleft) call drwcgm(xtip,ytip) xrite = xstart+del*cos(alpha+theta) yrite = ystart+del*sin(alpha+theta) call drwcgm(xrite,yrite) call drwcgm(xbase,ybase) return end subroutine arrgon ( arrow, xstart, ystart, xtip, ytip ) ! !*********************************************************************** ! !! ARRGON can be used to make a polygonal drawing of an arrow at any ! point on a graph. ! ! The arrow will stretch between two user specified points. ! ! The "head" of the arrow may be fatter or thinner than expected ! if the X and Y scales of the graph are not in the same ! proportions. ! ! It might be nice to include an "OUTLINE" option that draws the ! arrow outline in black, and fills it with the current color. ! ! By the way, the current color should either be a constant, ! or depend on the velocity magnitude, choosable by the user. ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! XSTART, ! YSTART Input, real XSTART, YSTART, the starting point for the ! arrow. ! ! XTIP, ! YTIP Input, real XTIP, YTIP, the end point for the arrow. ! ! ! left ! 5 ! |\ ! 7--------------6 \ ! start base 4tip ! 1--------------2 / ! |/ ! 3 ! rite ! real pi parameter (pi = 3.1415926) ! integer npts parameter (npts = 7) ! real alpha character ( len = 10 ) arrow real dist logical s_eqi real theta real xpts(npts+1) real xstart real xtip real ypts(npts+1) real ystart real ytip ! ! if ( xstart == xtip.and.ystart == ytip ) then return end if dist = sqrt((xtip-xstart)**2+(ytip-ystart)**2) theta = 0.5*pi-atan2(2.0,1.0) alpha = atan2(ytip-ystart,xtip-xstart) xpts(1) = xstart + dist*sin(alpha)/10.0 ypts(1) = ystart - dist*cos(alpha)/10.0 xpts(2) = xstart + dist*sin(alpha)/10.0 + dist*cos(alpha)*2.0/3.0 ypts(2) = ystart - dist*cos(alpha)/10.0 + dist*sin(alpha)*2.0/3.0 xpts(3) = xstart + dist*sqrt(5.0)*cos(alpha-theta)/3.0 ypts(3) = ystart + dist*sqrt(5.0)*sin(alpha-theta)/3.0 xpts(4) = xtip ypts(4) = ytip xpts(5) = xstart + dist*sqrt(5.0)*cos(alpha+theta)/3.0 ypts(5) = ystart + dist*sqrt(5.0)*sin(alpha+theta)/3.0 xpts(6) = xstart - dist*sin(alpha)/10.0 + dist*cos(alpha)*2.0/3.0 ypts(6) = ystart + dist*cos(alpha)/10.0 + dist*sin(alpha)*2.0/3.0 xpts(7) = xstart - dist*sin(alpha)/10.0 ypts(7) = ystart + dist*cos(alpha)/10.0 xpts(8) = xpts(1) ypts(8) = ypts(1) if ( s_eqi ( arrow,'hollow') ) then call plylin(npts+1,xpts,ypts) end if if ( s_eqi ( arrow,'solid') ) then call plygon(npts,xpts,ypts) end if return end subroutine bfl ( ielem, in, w, dwdx, dwdy, maxnpe, nelem, node, np, xc, & xq, yc, yq ) ! !*********************************************************************** ! !! BFL evaluates a particular linear basis function at a point ! in a nonisoparametric element. ! ! ^ ! | 2 ! | /| ! Y | / | ! | / | ! | 1---3 ! | ! +------------> ! X ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! IELEM Input, integer IELEM, the number of the element we are ! examining. This will be a value between 1 and NELEM. ! ! IN Input, integer IN, the number of the basis function we ! want. This will be a value between 1 and 3. ! ! W, ! DWDX, ! DWDY Output, real W, DWDX, DWDY, the value of the ! IN-th basis function and its X and Y derivatives, at the ! given point. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! Input, integer NELEM, the number of elements. ! ! NODE Input, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Input, integer NP, the number of nodes. ! ! XC Input, real XC(NP), the X coordinates of the ! nodes. ! ! XQ Input, real XQ, the X coordinate of the point ! where the basis function is to be evaluated. ! ! YC Input, real YC(NP), the Y coordinates of the nodes ! ! YQ Input, real YQ, the Y coordinate of the point wher ! the basis function is to be evaluated. ! integer maxnpe integer nelem integer np ! real d real dwdx real dwdy integer i1 integer i2 integer i3 integer ielem integer in integer in1 integer in2 integer in3 integer node(maxnpe,nelem) real w real xc(np) real xq real yc(np) real yq ! ! IN1, IN2, and IN3 are the local node numbers of the three ! corner nodes. ! in1 = in in2 = mod(in,3)+1 in3 = mod(in+1,3)+1 ! ! I1, I2, I3 are the global node numbers. ! i1 = node(in1,ielem) i2 = node(in2,ielem) i3 = node(in3,ielem) d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) if ( d == 0.0 ) then write(*,*)' ' write(*,*)'BFL - Fatal error!' write(*,*)' D = 0' write(*,*)' Element IELEM = ',ielem write(*,*)' I1, XC(I1), YC(I1) = ',i1,xc(i1),yc(i1) write(*,*)' I2, XC(I2), YC(I2) = ',i2,xc(i2),yc(i2) write(*,*)' I3, XC(I3), YC(I3) = ',i3,xc(i3),yc(i3) stop end if dwdx = (yc(i2)-yc(i3))/d dwdy = (xc(i3)-xc(i2))/d w = 1.0 + dwdx*(xq-xc(i1)) + dwdy*(yq-yc(i1)) return end subroutine bfq ( ielem,in,w,dwdx,dwdy,maxnpe,nelem,node,np,xc,xq,yc,yq ) ! !*********************************************************************** ! !! BFQ evaluates a particular quadratic basis function at a point ! in a nonisoparametric element. ! ! ^ ! | 2 ! | /| ! Y | 4 5 ! | / | ! | 1-6-3 ! | ! +------------> ! X ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! IELEM Input, integer IELEM, the number of the element we are ! examining. This will be a value between 1 and NELEM. ! ! IN Input, integer IN, the number of the basis function we ! want. This will be a value between 1 and 6. Functions ! 1 through 3 are associated with corners, 4 though 6 ! with sides. ! ! W, ! DWDX, ! DWDY Output, real W, DWDX, DWDY, the value of the ! IN-th basis function and its X and Y derivatives, at the ! given point. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! Input, integer NELEM, the number of elements. ! ! NODE Input, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Input, integer NP, the number of nodes. ! ! XC Input, real XC(NP), the X coordinates of the ! nodes. ! ! XQ Input, real XQ, the X coordinate of the point ! where the basis function is to be evaluated. ! ! YC Input, real YC(NP), the Y coordinates of the nodes ! ! YQ Input, real YQ, the Y coordinate of the point wher ! the basis function is to be evaluated. ! integer maxnpe integer nelem integer np ! real c real d real dwdx real dwdy integer i1 integer i2 integer i3 integer ielem integer in integer in1 integer in2 integer in3 integer node(maxnpe,nelem) real s real t real w real xc(np) real xq real yc(np) real yq ! ! Case 1: We are inquiring about a basis function associated ! with a corner. ! ! Notice that the basis function W is zero exactly if ! T is 0 or T is 1/2. ! ! IN1, IN2, and IN3 are the local node numbers of the three ! corner nodes, and I1, I2 and I3 are the corresponding ! global node numbers, which are used to look up the X and ! Y coordinates of the nodes. ! if ( 1 <= in.and.in <= 3 ) then in1 = in in2 = mod(in,3)+1 in3 = mod(in+1,3)+1 i1 = node(in1,ielem) i2 = node(in2,ielem) i3 = node(in3,ielem) d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) if ( d == 0.0 ) then write(*,*)' ' write(*,*)'BFQ - Fatal error!' write(*,*)' D = 0' write(*,*)' Element IELEM = ',ielem write(*,*)' I1, XC(I1), YC(I1) = ',i1,xc(i1),yc(i1) write(*,*)' I2, XC(I2), YC(I2) = ',i2,xc(i2),yc(i2) write(*,*)' I3, XC(I3), YC(I3) = ',i3,xc(i3),yc(i3) stop end if t = 1.0+( (xq -xc(i1))*(yc(i2)-yc(i3)) & +(xc(i3)-xc(i2))*(yq -yc(i1)) )/d w = t*(2.0*t-1.0) dwdx = (yc(i2)-yc(i3))*(4.0*t-1.0)/d dwdy = (xc(i3)-xc(i2))*(4.0*t-1.0)/d ! ! Case 2: We are inquiring about a basis function associated ! with a midpoint. ! else if ( in >= 4.and.in <= 6 ) then in1 = in-3 in2 = mod(in-3,3)+1 in3 = mod(in-2,3)+1 i1 = node(in1,ielem) i2 = node(in2,ielem) i3 = node(in3,ielem) d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1)) & -(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) if ( d == 0.0 ) then write(*,*)' ' write(*,*)'BFQ - Fatal error!' write(*,*)' D = 0' write(*,*)' Element IELEM = ',ielem write(*,*)' I1, XC(I1), YC(I1) = ',i1,xc(i1),yc(i1) write(*,*)' I2, XC(I2), YC(I2) = ',i2,xc(i2),yc(i2) write(*,*)' I3, XC(I3), YC(I3) = ',i3,xc(i3),yc(i3) stop end if c = (xc(i3)-xc(i2))*(yc(i1)-yc(i2))-(xc(i1)-xc(i2))*(yc(i3)-yc(i2)) if ( c == 0.0 ) then write(*,*)' ' write(*,*)'BFQ - Fatal error!' write(*,*)' C = 0' write(*,*)' Element IELEM = ',ielem write(*,*)' I1, XC(I1), YC(I1) = ',i1,xc(i1),yc(i1) write(*,*)' I2, XC(I2), YC(I2) = ',i2,xc(i2),yc(i2) write(*,*)' I3, XC(I3), YC(I3) = ',i3,xc(i3),yc(i3) stop end if t = 1.0+( (xq -xc(i1))*(yc(i2)-yc(i3))+(xc(i3)-xc(i2))*(yq -yc(i1)) )/d s = 1.0+( (xq -xc(i2))*(yc(i3)-yc(i1))+(xc(i1)-xc(i3))*(yq -yc(i2)) )/c w = 4.0 * s*t dwdx = 4.0 * ((yc(i3)-yc(i1))*t/c + (yc(i2)-yc(i3))*s/d) dwdy = 4.0 * ((xc(i1)-xc(i3))*t/c + (xc(i3)-xc(i2))*s/d) else write(*,*)' ' write(*,*)'BFQ - Fatal error!' write(*,*)' Request for basis function IN = ',in write(*,*)' but IN must be between 1 and 6.' stop end if return end subroutine bfrefl ( w,dwdx,dwdy,detadx,detady,dxsidx,dxsidy,eta,iq,xsi ) ! !*********************************************************************** ! !! BFREFL evaluates one of the three linear basis functions, ! and its X and Y derivatives, at a particular point (X,Y) ! in a particular element, by referring to the corresponding ! points (XSI,ETA) in the reference triangle. ! ! It is assumed that we already know the value of the jacobian ! of the isoparametric transformation between the (XSI, ETA) and ! (X, Y) spaces. The four entries of the jacobian are ! symbolically named DETADX, DETADY, DXSIDX and DXSIDY, and ! we know that the jacobian gives us the following relation ! between derivatives with respect to XSI and ETA, and derivatives ! with respect to X and Y: ! ! dF/dX = dF/dXsi dXsi/dX + dF/dEta dEta/dX ! dF/dY = dF/dXsi dXsi/dY + dF/dEta dEta/dY ! ! Here is a graph of the (XSI, ETA) reference triangle we will ! use. ! ! ^ ! | ! 1 + 2 ! | /| ! ETA | / | ! | / | ! 0 + 1---3 ! | ! +----+---+---> ! 0 1 ! ! XSI ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! W, ! DWDX, ! DWDY Output, real W, DWDX, DWDY, the value of the basis ! function, and its derivatives with respect to X and Y, at ! the point (ETA,XSI). ! ! DETADX, ! DETADY Input, real DETADX, DETADY, the partial derivative ! d ETA/d X and d ETA/d Y at (ETA,XSI). ! ! DXSIDX, ! DXSIDY Input, real DXSIDX, DXSIDY, the partial derivative ! d XSI/d X and d XSI/d Y at (ETA,XSI). ! ! ETA Input, real ETA, the local ETA coordinate ! at which the basis information is desired. ! ! IQ Input, integer IQ, the local node number, between 1 and ! 3, whose basis function is being evaluated. ! ! XSI Input, real XSI, the local XSI coordinate ! at which the basis information is desired. ! real detadx real detady real dwdeta real dwdx real dwdxsi real dwdy real dxsidx real dxsidy real eta integer iq real w real xsi ! if ( iq == 1 ) then w = 1.0-xsi dwdxsi = -1.0 dwdeta = 0.0 else if ( iq == 2 ) then w = eta dwdxsi = 0.0 dwdeta = 1.0 else if ( iq == 3 ) then w = xsi-eta dwdxsi = 1.0 dwdeta = -1.0 else write(*,*)' ' write(*,*)'BFRefL - Fatal error!' write(*,*)' Request for basis function IQ = ',iq write(*,*)' but IQ must be between 1 and 3.' stop end if dwdx = dwdxsi*dxsidx+dwdeta*detadx dwdy = dwdxsi*dxsidy+dwdeta*detady return end subroutine bfrefq ( w,dwdx,dwdy,detadx,detady,dxsidx,dxsidy,eta,iq,xsi ) ! !*********************************************************************** ! !! BFREFQ evaluates one of the six quadratic basis functions, ! and its X and Y derivatives, at a particular point in a ! finite element, by referring to the reference triangle. ! ! The point we are interested in is referred to by its coordinates ! in the reference triangle. That is, we are given coordinates ! (XSI, ETA), even though, physically, we are interested ! in points in (X, Y) space. ! ! It is assumed that we already know the value of the jacobian ! of the isoparametric transformation between the (XSI, ETA) and ! (X, Y) spaces. The four entries of the jacobian are ! symbolically named DETADX, DETADY, DXSIDX and DXSIDY, and ! we know that the jacobian gives us the following relation ! between derivatives with respect to XSI and ETA, and derivatives ! with respect to X and Y: ! ! d F(X,Y)/dX (d XSI/dX d ETA/dX ) ( d F(XSI, ETA)/d XSI ) ! d F(X,Y)/dY = (d XSI/dY d ETA/dY ) * ( d F(XSI, ETA)/d ETA ) ! ! Here is a graph of the (XSI, ETA) reference triangle we will ! use. ! ! ^ ! | ! 1 + 2 ! | /| ! ETA | 4 5 ! | / | ! 0 + 1-6-3 ! | ! +----+---+---> ! 0 1 ! ! XSI ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! W, ! DWDX, ! DWDY Output, real W, DWDX, DWDY, the value of the basis ! function, and its derivatives with respect to X and Y, at ! the point (XSI,ETA). ! ! DETADX, ! DETADY Input, real DETADX, DETADY, the partial derivatives ! d ETA/d X and d ETA/d Y at (XSI,ETA). ! ! DXSIDX, ! DXSIDY Input, real DXSIDX, DXSIDY, the partial derivatives ! d XSI/d X and d XSI/d Y at (XSI,ETA). ! ! ETA Input, real ETA, the ETA coordinate of the point. ! ! IQ Input, integer IQ, the local node number, between 1 and 6, ! whose basis function is being evaluated. ! ! XSI Input, real XSI, the XSI coordinate of the point. ! real detadx real detady real dwdeta real dwdx real dwdxsi real dwdy real dxsidx real dxsidy real eta integer iq real w real xsi ! ! Find the formula for the desired basis function, and evaluate it. ! Also evaluate the partial derivatives d/d XSI and d/d ETA. ! ! Basis function 1 is zero if XSI = 0.5 or XSI=1. ! if ( iq == 1 ) then w = (2.0*xsi-1.0) * (xsi-1.0) dwdxsi = -3.0+4.0*xsi dwdeta = 0.0 ! ! Basis function 2 is zero if ETA = 0 or ETA=0.5. ! else if ( iq == 2 ) then w = eta * (2.0*eta-1.0) dwdxsi = 0.0 dwdeta = -1.0+4.0*eta ! ! Basis function 3 is zero if XSI = ETA, or XSI=ETA+0.5 ! else if ( iq == 3 ) then w = (xsi-eta) * (2.0*xsi-2.0*eta-1.0) dwdxsi = -1.0+4.0*xsi-4.0*eta dwdeta = 1.0-4.0*xsi+4.0*eta ! ! Basis function 4 is zero if ETA = 0 or XSI=1. ! else if ( iq == 4 ) then w = 4.0 * eta * (1.0-xsi) dwdxsi = -4.0*eta dwdeta = 4.0-4.0*xsi ! ! Basis function 5 is zero if ETA = 0 or XSI=ETA. ! else if ( iq == 5 ) then w = 4.0 * eta * (xsi-eta) dwdxsi = 4.0*eta dwdeta = 4.0*xsi-8.0*eta ! ! Basis function 6 is zero if XSI = ETA or XSI=1. ! else if ( iq == 6 ) then w = 4.0 * (xsi-eta) * (1.0-xsi) dwdxsi = 4.0-8.0*xsi+4.0*eta dwdeta = -4.0+4.0*xsi else write(*,*)' ' write(*,*)'BFRefQ - Fatal error!' write(*,*)' Illegal value of IQ = ',iq stop end if ! ! Convert the d/d XSI and d/d ETA derivatives to d/d X and d/d Y. ! dwdx = dwdxsi*dxsidx+dwdeta*detadx dwdy = dwdxsi*dxsidy+dwdeta*detady return end subroutine boundary_set ( eqn,jbound,maxbou,maxnpe,nbound,nelem,node,np,npe ) ! !*********************************************************************** ! !! BOUNDARY_SET deduces which edges of the finite element triangle constitute ! boundaries of the region. ! ! It does this stupidly, by considering each element, and then each ! edge of that element, and seeing if it occurs (in the opposite ! orientation) in any other element. ! ! For large regions, SETBND is bound to be slow. If you know ! where your boundaries occur, then you might want to replace ! BOUND by a simpler routine. However, SETBND is guaranteed ! to work for a wide variety of regions. ! ! If an edge only occurs in one element, and no velocities ! are specified along it, then it is a wall, and we record it ! as a boundary. ! ! Only edges formed by the three corner nodes are considered. ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! EQN Input, character ( len = 2 ) EQN(3,NP). ! EQN records the "type" of each equation that will be generated, and ! which is associated with an unknown. Note that most boundary ! conditions do not result in an equation. The current values are: ! ! 'U' The horizontal momentum equation. ! 'UW' The condition U = 0 applied at a node on a fixed wall. ! 'V' The vertical momentum equation. ! 'VW' The condition V = 0 applied at a node on a fixed wall. ! 'P' The continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! ! JBOUND Output, integer JBOUND(5,MAXBOU) ! ! For each line segment of the boundary: ! ! JBOUND(1,I) contains the element number; ! ! JBOUND(2,I) contains the local node number of one corner ! of the element, which forms the edge; ! ! JBOUND(2,I) contains the "next" node along the edge. ! If the element is linear, this is the other corner node. ! If the element is quadratic, this is the midside node along ! the edge. ! ! JBOUND(4,I) contains the "next" node along the edge. ! If the element is linear, this is 0. ! If the element is quadratic, this is the other corner node ! along the edge. ! ! JBOUND(5,I) contains: ! 0 if the boundary is a wall (U = V=0); ! 1 if the boundary is open. ! ! MAXBOU Input, integer MAXBOU. ! The amount of storage available for the IBOUND array. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! NBOUND Input, integer NBOUND. ! The number of points (XBOUND(I),YBOUND(I)) used to ! define the boundary. ! ! Input, integer NELEM. ! The number of elements. ! ! NODE Input, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Input, integer NP, the number of nodes. ! ! NPE Input, integer NPE. ! NPE is the number of nodes per element, which should be ! 3 for linear elements and 6 for quadratics. integer maxbou integer maxnpe integer nelem integer np ! character ( len = 2 ) eqn(3,np) integer i integer j integer jbound(5,maxbou) integer jp1 integer k integer l integer lp1 integer m1 integer m2 integer n1 integer n2 integer nbound integer nedge integer node(maxnpe,nelem) integer npe ! if ( npe == 3 ) then nedge = 3 else if ( npe == 4 ) then nedge = 4 else if ( npe == 6 ) then nedge = 3 end if nbound = 0 ! ! For each element I... ! do i = 1,nelem ! ! ...consider, for J = 1 to NEDGE, edges of the form ! (N1, N2) = (NODE(J,I), NODE(J+1,I))... ! do j = 1,nedge n1 = node(j,i) jp1 = j+1 if ( j == nedge)jp1 = 1 n2 = node(jp1,i) ! ! ...and now, for each element K,... ! do k = 1,nelem ! ! ...where K is different from I, ... ! if ( k /= i ) then ! ! ...consider, for L = 1 to NEDGE, edges of the form ! (M1, M2) = (NODE(L,K), NODE(L+1,K))... ! do l = 1,nedge m1 = node(l,k) lp1 = l+1 if ( l == nedge)lp1 = 1 m2 = node(lp1,k) ! ! And if (N1, N2) = (M1, M2), or (N1, N2) = (M2, M1), then skip out. ! if ( (m1 == n1.and.m2 == n2) .or. (m1 == n2.and.m2 == n1) ) then go to 10 end if end do end if end do ! ! Otherwise, the pair (N1, N2) occurs in no other element. Therefore, ! it is a boundary edge. ! if ( nbound >= maxbou ) then write(*,*)' ' write(*,*)'SETBND - Warning!' write(*,*)' Ran out of space to store boundary!' write(*,*)' Currently using ',nbound,' boundary edges.' return end if ! ! Increment the number of boundary edges found, and store ! the indices of the element, the first node, the middle node ! (if quadratic), and the second node. ! nbound = nbound+1 jbound(1,nbound) = i jbound(2,nbound) = j if ( npe == 3 ) then jbound(3,nbound) = jp1 else if ( npe == 4 ) then jbound(3,nbound) = jp1 else if ( npe == 6 ) then jbound(3,nbound) = j+3 end if if ( npe == 3 ) then jbound(4,nbound) = 0 else if ( npe == 4 ) then jbound(4,nbound) = 0 else if ( npe == 6 ) then jbound(4,nbound) = jp1 end if ! ! Now try to figure out what kind of boundary segment it is. ! ! It is OPEN if there is an unknown U velocity at either end node ! which does not correspond to a wall, ! or there is a specified U velocity at either end node. ! if ( eqn(1,n1) == 'UW'.and.eqn(1,n2) == 'UW' ) then jbound(5,nbound) = 0 else jbound(5,nbound) = 1 end if 10 continue end do end do write(*,*)' ' write(*,*)'SETBND - Note:' write(*,*)' Found NBOUND = ',nbound,' boundary edges.' return end subroutine boundary_show ( eflag, etaref, isotri, jbound, line, maxbou, maxnpe, & maxobj, nbound, nelem, nflag, node, np, npe, xc, xsiref, yc ) ! !*********************************************************************** ! !! BOUNDARY_SHOW displays the boundary. ! ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical EFLAG(NELEM), flags the elements to be displayed. ! If EFLAG(I) is TRUE, then element I should be displayed. ! ! ETAREF Input, real ETAREF(MAXNPE). ! The ETA coordinates of the nodes of the reference triangle. ! ! ISOTRI Input, integer ISOTRI(NELEM). ! ! 0, if element I is not isoparametric. ! 1 or 2, if element I is isoparametric. ! ! JBOUND Input, integer JBOUND(5,MAXBOU) ! ! For each line segment of the boundary: ! ! JBOUND(1,I) contains the element number; ! ! JBOUND(2,I) contains the local node number of one corner ! of the element, which forms the edge; ! ! JBOUND(2,I) contains the "next" node along the edge. ! If the element is linear, this is the other corner node. ! If the element is quadratic, this is the midside node along ! the edge. ! ! JBOUND(4,I) contains the "next" node along the edge. ! If the element is linear, this is 0. ! If the element is quadratic, this is the other corner node ! along the edge. ! ! JBOUND(5,I) contains: ! 0 if the boundary is a wall (U = V=0); ! 1 if the boundary is open. ! ! MAXBOU Input, integer MAXBOU. ! The amount of storage available for the IBOUND array. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! NBOUND Input, integer NBOUND. ! The number of points (XBOUND(I),YBOUND(I)) used to ! define the boundary. ! ! Input, integer NELEM, the number of elements. ! ! NFLAG Input, logical NFLAG(NP). ! ! NFLAG is used to "flag" which nodes are active, ! that is, to be displayed, and which not, in the graph. ! ! NODE Input, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Input, integer NP, the number of nodes. ! ! NPE Input, integer NPE. ! NPE is the number of nodes per element, which should be ! 3 for linear elements and 6 for quadratics. ! ! XC Input, real XC(NP). ! The X coordinates of the nodes. ! ! XSIREF Input, real XSIREF(MAXNPE). ! The XSI coordinates of the nodes of the reference ! triangle. ! ! YC Input, real YC(NP). ! The Y coordinates of the nodes. ! integer maxbou integer maxnpe integer maxobj integer nelem integer np ! logical eflag(nelem) real etaref(maxnpe) integer i integer ielem integer isotri(nelem) integer j integer j1 integer j2 integer jbound(5,maxbou) integer jhi integer line(maxobj) integer nbound logical nflag(np) integer ng1 integer ng2 integer ng3 integer node(maxnpe,nelem) integer npe real x1 real x2 real x3 real xc(np) real xsiref(maxnpe) real xtemp real y1 real y2 real y3 real yc(np) real ytemp ! ! Draw the line segments: ! for linears: from node N1 to node N2; ! for quadratics: from node N1 to N2, and from N2 to N3. ! ! ! Consider the I-th triangle edge that lies on the boundary. ! do i = 1,nbound ! ! Retrieve the element number IELEM of the triangle that contains ! the edge. ! ielem = jbound(1,i) ! ! If the element is visible, then proceed, else skip to next ! triangle edge. ! if ( eflag(ielem) ) then ! ! Boundary stretch is part of nonisoparametric element. ! if ( npe == 3 .or. npe == 4 .or. isotri(ielem) == 0 ) then ng1 = node(jbound(2,i),ielem) ng2 = node(jbound(3,i),ielem) if ( npe == 3 ) then ng3 = 0 else if ( npe == 4 ) then ng3 = 0 else ng3 = node(jbound(4,i),ielem) end if ! ! Draw a solid line, if requested. ! if ( line(1) == 0 .or. line(1) == 2 ) then if ( nflag(ng1) .and. nflag(ng2) ) then call movcgm(xc(ng1),yc(ng1)) call drwcgm(xc(ng2),yc(ng2)) end if if ( npe == 6 ) then if ( nflag(ng2).and.nflag(ng3) ) then call movcgm(xc(ng2),yc(ng2)) call drwcgm(xc(ng3),yc(ng3)) end if end if ! ! Draw a solid wall. ! else if ( jbound(5,i) == 0 ) then if ( nflag(ng1).and.nflag(ng2) ) then call movcgm(xc(ng1),yc(ng1)) call drwcgm(xc(ng2),yc(ng2)) end if if ( npe == 6 ) then if ( nflag(ng2).and.nflag(ng3) ) then call movcgm(xc(ng2),yc(ng2)) call drwcgm(xc(ng3),yc(ng3)) end if end if ! ! Draw an dashed open boundary. ! else if ( nflag(ng1).and.nflag(ng2) ) then call movcgm(xc(ng1),yc(ng1)) xtemp = xc(ng1)+0.25*(xc(ng2)-xc(ng1)) ytemp = yc(ng1)+0.25*(yc(ng2)-yc(ng1)) call drwcgm(xtemp,ytemp) xtemp = xc(ng1)+0.75*(xc(ng2)-xc(ng1)) ytemp = yc(ng1)+0.75*(yc(ng2)-yc(ng1)) call movcgm(xtemp,ytemp) call drwcgm(xc(ng2),yc(ng2)) end if if ( npe == 6 ) then if ( nflag(ng2).and.nflag(ng3) ) then call movcgm(xc(ng2),yc(ng2)) xtemp = xc(ng2)+0.25*(xc(ng3)-xc(ng2)) ytemp = yc(ng2)+0.25*(yc(ng3)-yc(ng2)) call drwcgm(xtemp,ytemp) xtemp = xc(ng2)+0.75*(xc(ng3)-xc(ng2)) ytemp = yc(ng2)+0.75*(yc(ng3)-yc(ng2)) call movcgm(xtemp,ytemp) call drwcgm(xc(ng3),yc(ng3)) end if end if end if end if ! ! Boundary stretch that is part of an isoparametric element. ! else if ( npe == 6 ) then jhi = 2 end if do j = 1,jhi j1 = jbound(j+1,i) j2 = jbound(j+2,i) x1 = xsiref(j1) y1 = etaref(j1) x2 = xsiref(j2) y2 = etaref(j2) if ( line(1) == 0 .or. line(1) == 2 ) then call isoln6(x1,y1,x2,y2,ielem,maxnpe,nelem,node,np,npe,xc,yc) else if ( jbound(5,i) == 0 ) then call isoln6(x1,y1,x2,y2,ielem,maxnpe,nelem,node,np,npe,xc,yc) else x1 = xsiref(j1) y1 = etaref(j1) x3 = x1+0.25*(x2-x1) y3 = y1+0.25*(y2-y1) call isoln6(x1,y1,x3,y3,ielem,maxnpe,nelem,node,np,npe,xc,yc) x3 = x1+0.75*(x2-x1) y3 = y1+0.75*(y2-y1) call isoln6(x3,y3,x2,y2,ielem,maxnpe,nelem,node,np,npe,xc,yc) end if end if end do end if end if end do return end subroutine box ( xmin, xmax, ymin, ymax ) ! !******************************************************************************* ! !! BOX draws a rectangle whose corners are specified by the user. ! ! ! Discussion: ! ! The rectangle drawn by box has the corners: ! ! (XMIN,YMAX) (XMAX,YMAX) ! ! (XMIN,YMIN) (XMAX,YMIN) ! ! BOX can be used to place a rectangle anywhere in the picture. ! However, BOX may also be used to place a rectangle around the ! entire picture, producing a "frame". ! ! The DRAWCGM routine PLYLIN is used to draw the box, and hence ! the color of the line may be set by calling the DRAWCGM routine ! LINCLR first. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real XMIN, XMAX, the minimum and maximum X ! coordinates of the box. ! ! Input, real YMIN, YMAX, the minimum and maximum Y ! coordinates of the box. ! integer, parameter :: npoints = 5 ! real x(npoints) real xmax real xmin real y(npoints) real ymax real ymin ! x(1) = xmin y(1) = ymin x(2) = xmax y(2) = ymin x(3) = xmax y(3) = ymax x(4) = xmin y(4) = ymax x(5) = xmin y(5) = ymin call plylin ( npoints, x, y ) return end subroutine buzz ( dev, x1min, x1max, y1min, y1max ) ! !******************************************************************************* ! !! BUZZ is just "busy work" to force XWS to draw the whole picture. ! ! ! Discussion: ! ! BUZZ seems to fix a problem that occurs when the XWS interface is ! used. In that case, the last bit of the graph is not drawn. So ! here, we just make the last bit of the graph something we don't ! care about. ! ! Modified: ! ! 01 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) DEV, the graphics output device: ! cgmb - CGM binary file. ! ps - PostScript file. ! XWS - X window screen (interactive). ! ! Input, real X1MIN, X1MAX, the minimum and maximum X ! coordinates of the plot, which includes a small grace margin. ! ! Input, real Y1MIN, Y1MAX, the minimum and maximum Y ! coordinates of the plot, which includes a small grace margin. ! character ( len = * ) dev integer i integer icolor logical s_eqi real x1max real x1min real y1max real y1min ! if ( s_eqi ( dev, 'XWS' ) ) then icolor = 0 icolor = 1 call linclr ( icolor ) do i = 1, 10000 call box ( x1min, x1max, y1min, y1max ) end do end if return end subroutine ch_cap ( c ) ! !******************************************************************************* ! !! CH_CAP capitalizes a single character. ! ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! character c integer itemp ! itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end subroutine cbar ( icolor,jcmax,jcmin,maxobj,ncon,smax,smin,srange,x1,x2,y1,y2) ! !*********************************************************************** ! !! CBAR draws a color bar in the rectangle whose corners are ! (X1,Y1) and (X2,Y2). ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! ICOLOR Input, integer ICOLOR(MAXOBJ). ! ICOLOR contains the color index for each object. ! ! JCMAX, ! JCMIN Input, integer JCMAX, JCMIN, the maximum and ! minimum color indices to use in the color bar. ! ! MAXOBJ Input, integer MAXOBJ. ! The number of graphical "objects". ! ! NCON Input, integer NCON, the number of color contour ! regions drawn, and hence, the number of colors ! to be displayed in the color bar. ! ! SMAX, ! SMIN Input, real SMAX, SMIN, the maximum and minimum ! values of the quantity whose color contours are ! being drawn. These numbers will be printed along ! with the color bar. ! ! SRANGE Input, real SRANGE. ! The maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! X1, ! X2, ! Y1, ! Y2 Input, real X1, X2, Y1, Y2, specify the minimum and ! maximum X and Y coordinates of the color bar. ! integer maxobj ! real angle character ( len = 14 ) chrrel character ( len = 14 ) chrtmp real cwide character ( len = 6 ) flush integer i integer icolor(maxobj) integer jcmax integer jcmin integer jcolor integer lent integer ncon real pwide real smax real smin real srange real x real x1 real x2 real xcorn(5) real xl real xr real y real y1 real y2 real ycorn(5) ! ! ycorn(1) = y1 ycorn(2) = y1 ycorn(3) = y2 ycorn(4) = y2 ycorn(5) = y1 call linclr(icolor(1)) do i = 0,ncon xl = ((ncon+1-i)*x1+i*x2)/real(ncon+1) xr = ((ncon-i)*x1+(i+1)*x2)/real(ncon+1) xcorn(1) = xl xcorn(2) = xr xcorn(3) = xr xcorn(4) = xl xcorn(5) = xl jcolor = int(((ncon-i)*jcmin+i*jcmax)/real(ncon)) call filclr(jcolor) call plygon(4,xcorn,ycorn) call plylin(5,xcorn,ycorn) end do ! ! Print labels for the lowest and highest contours. ! cwide = 0.9*srange/40.0 chrtmp = chrrel(smin) call s_blank_delete ( chrtmp) lent = len_trim ( chrtmp ) angle = 0.0 x = x1 y = y1-1.5*cwide pwide = srange flush = 'left' call s_plot(angle,cwide,pwide,chrtmp(1:lent),x,y,flush) chrtmp = chrrel(smax) call s_blank_delete ( chrtmp) lent = len_trim ( chrtmp ) angle = 0.0 x = x2 y = y1-1.5*cwide pwide = srange flush = 'right' call s_plot(angle,cwide,pwide,chrtmp(1:lent),x,y,flush) return end subroutine cbox ( grace ) ! !*********************************************************************** ! !! CBOX draws a 16 by 16 box, filling each entry with a color ! from the current color table. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real GRACE. ! The size of the "grace" margin on the plot. ! integer npoly parameter (npoly = 5) ! real grace integer i integer icolor integer ierror integer j real x(npoly) real x1max real x1min real x2max real x2min real y(npoly) real y1max real y1min real y2max real y2min ! ! Set the coordinate system to be 0 < = X <= 16.0, ! and 0.0 < = Y <= 16.0 ! x2min = 0.0 x2max = 16.0 y2min = 0.0 y2max = 16.0 x1min = x2min-grace*(x2max-x2min) x1max = x2max+grace*(x2max-x2min) y1min = y2min-grace*(y2max-y2min) y1max = y2max+grace*(y2max-y2min) call setwcd(x1min,y1min,x1max,y1max,ierror) ! ! Draw the color boxes. ! icolor = 0 do i = 1,16 y(1) = 16-i y(2) = 16-i y(3) = 17-i y(4) = 17-i do j = 1,16 call filclr(icolor) icolor = icolor+1 x(1) = j-1 x(2) = j x(3) = j x(4) = j-1 call plygon(4,x,y) end do end do ! ! Draw black lines around the boxes. ! icolor = 1 call linclr(icolor) do i = 1,16 y(1) = 16-i y(2) = 17-i y(3) = 17-i y(4) = 16-i y(5) = y(1) do j = 1,16 x(1) = j-1 x(2) = j-1 x(3) = j x(4) = j x(5) = x(1) call plylin(5,x,y) end do end do return end function chrint ( intval ) ! !*********************************************************************** ! !! CHRINT accepts an integer and returns in CHRINT the 6-character ! representation of the integer, right justified, or '******' if ! the integer is too large or negative to fit in six positions. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! INTVAL Input, integer INTVAL, an integer variable to be ! converted. ! ! CHRINT Output (through function value), character*6 CHRINT, a 6 ! character representation of the integer, right justified. ! Thus, if INTVAL = 1,CHRINT=' 1'. CHRINT must be ! declared "character CHRINT*6" in the calling program. ! character ( len = 6 ) chrint character ( len = 6 ) chrtmp integer intval ! if ( intval > 999999 ) then chrtmp = '******' else if ( intval < -99999 ) then chrtmp = '-*****' else write(chrtmp,'(i6)')intval end if chrint = chrtmp return end function chrrel ( rval ) ! !*********************************************************************** ! !! CHRREL accepts a real number in RVAL and returns in CHRREL a ! 14-character right-justified representation of that number. ! ! ! Modified: ! ! 03 April 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! RVAL Input, real RVAL, a real number. ! ! CHRREL Output (through function value), character ( len = 14 ) CHRREL, ! a right-justified character variable containing the ! representation of RVAL, using a G14.6 format. ! character ( len = 14 ) chrrel character ( len = 14 ) chrtmp real rval ! ! ! We can't seem to write directly into CHRREL because of compiler ! quibbles. ! if ( real(int(rval)) == rval.and.abs(rval) < 1.0E+13 ) then write(chrtmp,'(i14)')int(rval) else write(chrtmp,'(g14.6)')rval end if chrrel = chrtmp return end subroutine conc3 ( eflag,jcmax,jcmin,maxnpe,ncon,nelem,node,np,s, & smax,smin,xc,yc) ! !*********************************************************************** ! !! CONC3 uses color to indicate all the points which have a function ! value between given limits. ! ! CONC3 is used for quantities associated with a 3 node finite element. ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! EFLAG Input, logical EFLAG(NELEM). ! EFLAG is used to "flag" which elements are to be displayed. ! If EFLAG(I) is TRUE, then element I should be displayed. ! ! ETAREF Input, real ETAREF(MAXNPE). ! The ETA coordinates of the nodes of the reference ! triangle. ! ! ISOTRI Output, integer ISOTRI(NELEM). ! ! 0, if element I is not isoparametric. ! 1 or 2, if element I is isoparametric. ! ! JCMAX, ! JCMIN Input, integer JCMAX, JCMIN, the maximum and minimum color ! indices to use for contours. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! NCON Input, integer NCON, the number of color contours to use. ! ! Input, integer NELEM. ! The number of elements. ! ! NODE Input, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Input, integer NP, the number of nodes. ! ! NPE Input, integer NPE. ! NPE is the number of nodes per element, which should be ! 3 for linear elements and 6 for quadratics. ! ! S Input, real S(NP). ! S is a scalar quantity associated with the nodes. ! This routine only looks at the values associated with ! corner element nodes. ! ! SMAX, ! SMIN Input, real SMAX, SMIN, the maximum and minimum values of S. ! ! XC Input, real XC(NP). ! The X coordinates of the nodes. ! ! XSIREF Input, real XSIREF(MAXNPE). ! The XSI coordinates of the nodes of the reference ! triangle. ! ! YC Input, real YC(NP). ! The Y coordinates of the nodes. ! integer maxnpe integer nelem integer np ! logical eflag(nelem) integer i integer i1 integer i2 integer i3 integer jcmax integer jcmin integer ncon integer node(maxnpe,nelem) real s(np) real smax real smin real xc(np) real yc(np) ! do i = 1,nelem if ( eflag(i) ) then i1 = node(1,i) i2 = node(2,i) i3 = node(3,i) call tricno(i1,i2,i3,jcmax,jcmin,ncon,np,s,smax,smin,xc,yc) end if end do return end subroutine conc4 ( eflag,jcmax,jcmin,maxnpe,ncon,nelem,node,np,s, & smax,smin,xc,yc) ! !*********************************************************************** ! !! CONC4 uses color to indicate all the points which have a function ! value between given limits. ! ! CONC4 is used for quantities associated with a nonisoparametric ! rectangular 4 node finite element. ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! EFLAG Input, logical EFLAG(NELEM). ! EFLAG is used to "flag" which elements are to be displayed. ! If EFLAG(I) is TRUE, then element I should be displayed. ! ! ETAREF Input, real ETAREF(MAXNPE). ! The ETA coordinates of the nodes of the reference ! triangle. ! ! ISOTRI Output, integer ISOTRI(NELEM). ! ! 0, if element I is not isoparametric. ! 1 or 2, if element I is isoparametric. ! ! JCMAX, ! JCMIN Input, integer JCMAX, JCMIN, the maximum and minimum color ! indices to use for contours. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! NCON Input, integer NCON, the number of color contours to use. ! ! Input, integer NELEM. ! The number of elements. ! ! NODE Input, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Input, integer NP, the number of nodes. ! ! NPE Input, integer NPE. ! NPE is the number of nodes per element, which should be ! 3 for linear elements and 6 for quadratics. ! ! S Input, real S(NP). ! S is a scalar quantity associated with the nodes. ! This routine only looks at the values associated with ! corner element nodes. ! ! SMAX, ! SMIN Input, real SMAX, SMIN, the maximum and minimum values of S. ! ! XC Input, real XC(NP). ! The X coordinates of the nodes. ! ! XSIREF Input, real XSIREF(MAXNPE). ! The XSI coordinates of the nodes of the reference ! triangle. ! ! YC Input, real YC(NP). ! The Y coordinates of the nodes. ! integer maxnpe integer nelem integer np ! logical eflag(nelem) integer i integer i1 integer i2 integer i3 integer jcmax integer jcmin integer ncon integer node(maxnpe,nelem) real s(np) real smax real smin real xc(np) real yc(np) ! do i = 1,nelem if ( eflag(i) ) then i1 = node(1,i) i2 = node(2,i) i3 = node(3,i) call tricno(i1,i2,i3,jcmax,jcmin,ncon,np,s,smax,smin,xc,yc) i1 = node(3,i) i2 = node(4,i) i3 = node(1,i) call tricno(i1,i2,i3,jcmax,jcmin,ncon,np,s,smax,smin,xc,yc) end if end do return end subroutine conc6 ( eflag,etaref,isotri,jcmax,jcmin,maxnpe,ncon, & nelem,node,np,npe,s,smax,smin,xc,xsiref,yc ) ! !*********************************************************************** ! !! CONC6 uses color to indicate all the points which have a ! function value greater than a given value. ! ! CONC6 is used for quantities associated with the six ! corner nodes of a 6 node finite element, in particular, ! vorticity or velocity magnitude. ! ! CONC6 breaks this triangular element up into four three ! node triangles for further treatment. ! ! ! 2 ! /| ! / | ! 4--5 ! /| /| ! / |/ | ! 1--6--3 ! ! Thus, the four triangles are defined by the nodes as follows: ! ! Triangle Nodes ! ! 1 1, 4, 6 ! 2 2, 5, 4 ! 3 3, 6, 5 ! 4 4, 5, 6 ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! EFLAG Input, logical EFLAG(NELEM). ! EFLAG is used to "flag" which elements are to be displayed. ! If EFLAG(I) is TRUE, then element I should be displayed. ! ! ETAREF Input, real ETAREF(MAXNPE). ! The ETA coordinates of the nodes of the reference ! triangle. ! ! ISOTRI Output, integer ISOTRI(NELEM). ! ! 0, if element I is not isoparametric. ! 1 or 2, if element I is isoparametric. ! ! JCMAX, ! JCMIN Input, integer JCMAX, JCMIN, the maximum and minimum color ! indices to use for contours. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! NCON Input, integer NCON, the number of contours to draw. ! ! Input, integer NELEM. ! The number of elements. ! ! NODE Input, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Input, integer NP, the number of nodes. ! ! NPE Input, integer NPE. ! NPE is the number of nodes per element, which should be ! 3 for linear elements and 6 for quadratics. ! ! S Input, real S(NP), the value of S at each node. ! ! SMAX, ! SMIN Input, real SMAX, SMIN, the maximum and minimum values of S. ! ! XC Input, real XC(NP), the X coordinates of the nodes. ! ! XSIREF Input, real XSIREF(MAXNPE). ! The XSI coordinates of the nodes of the reference ! triangle. ! ! YC Input, real YC(NP), the Y coordinates of the nodes. ! integer maxnpe integer nelem integer np ! logical eflag(nelem) real eta1 real eta2 real eta3 real etaref(maxnpe) integer i integer i1 integer i2 integer i3 integer ielem integer isotri(nelem) integer j integer jcmax integer jcmin integer ncon integer node(maxnpe,nelem) integer npe integer nsub(4,3) real s(np) real s1 real s2 real s3 real smax real smin real xc(np) real xsi1 real xsi2 real xsi3 real xsiref(maxnpe) real yc(np) ! nsub(1,1) = 1 nsub(1,2) = 4 nsub(1,3) = 6 nsub(2,1) = 2 nsub(2,2) = 5 nsub(2,3) = 4 nsub(3,1) = 3 nsub(3,2) = 6 nsub(3,3) = 5 nsub(4,1) = 4 nsub(4,2) = 5 nsub(4,3) = 6 ! ! Draw the contour line by searching over each element. ! do i = 1,nelem ielem = i if ( eflag(i) ) then if ( npe == 3.or.isotri(i) == 0 ) then do j = 1,4 i1 = node(nsub(j,1),i) i2 = node(nsub(j,2),i) i3 = node(nsub(j,3),i) call tricno(i1,i2,i3,jcmax,jcmin,ncon,np,s,smax,smin,xc,yc) end do else do j = 1,4 i1 = node(nsub(j,1),i) i2 = node(nsub(j,2),i) i3 = node(nsub(j,3),i) eta1 = etaref(nsub(j,1)) eta2 = etaref(nsub(j,2)) eta3 = etaref(nsub(j,3)) s1 = s(i1) s2 = s(i2) s3 = s(i3) xsi1 = xsiref(nsub(j,1)) xsi2 = xsiref(nsub(j,2)) xsi3 = xsiref(nsub(j,3)) call tricis(eta1,eta2,eta3,ielem,jcmax,jcmin,maxnpe, & ncon,nelem,node,np,npe,s1,s2,s3,smax,smin,xc,xsi1, & xsi2,xsi3,yc) end do end if end if end do return end subroutine conc63 ( eflag,etaref,isotri,jcmax,jcmin,maxnpe,ncon, & nelem,node,np,npe,s,smax,smin,xc,xsiref,yc) ! !*********************************************************************** ! !! CONC63 uses color to indicate all the points which have a function ! value between given limits. ! ! CONC63 is used for quantities associated with the three corner nodes ! of a 6 node finite element, in particular, pressure. ! ! Modified: ! ! 12 January 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! EFLAG Input, logical EFLAG(NELEM). ! EFLAG is used to "flag" which elements are to be displayed. ! If EFLAG(I) is TRUE, then element I should be displayed. ! ! ETAREF Input, real ETAREF(MAXNPE). ! The ETA coordinates of the nodes of the reference ! triangle. ! ! ISOTRI Output, integer ISOTRI(NELEM). ! ! 0, if element I is not isoparametric. ! 1 or 2, if element I is isoparametric. ! ! JCMAX, ! JCMIN Input, integer JCMAX, JCMIN, the maximum and minimum color ! indices to use for contours. ! ! MAXNPE Input, integer MAXNPE. ! MAXNPE is the maximum number of nodes per element. ! ! NCON Input, integer NCON, the number of color contours to use. ! ! Input, integer NELEM. ! The number of elements. ! ! NODE Input, integer NODE(MAXNPE,MAXELM), or NODE(MAXNPE,NELEM). ! NODE contains, for each element, the global node numbers ! of its NPE nodes. ! ! For linear elements (NPE = 3), the order of the nodes probably ! doesn't matter, but we will draw them this way: ! ! 2 ! /| ! / | ! / | ! / | ! 1-------3 ! ! For quadratic elements (NPE = 6), the nodes must be given in a ! particular order, which is as follows: ! ! 2 ! /| ! / | ! 4 5 ! / | ! 1---6---3 ! ! Input, integer NP, the number of nodes. ! ! NPE Input, integer NPE. ! NPE is the number of nodes per element, which should be ! 3 for linear elements and 6 for quadratics. ! ! S Input, real S(NP). ! S is a scalar quantity associated with the nodes. ! This routine only looks at the values associated with ! corner element nodes. ! ! SMAX, ! SMIN Input, real SMAX, SMIN, the maximum and minimum values of S. ! ! XC Input, real XC(NP). ! The X coordinates of the nodes. ! ! XSIREF Input, real XSIREF(MAXNPE). ! The XSI coordinates of the nodes of the reference ! triangle. ! ! YC Input, real YC(NP). ! The Y coordinates of the nodes.