program main !*****************************************************************************80 ! !! MAIN is the main program for DISPLAY4. ! ! Discussion: ! ! DISPLAY4 displays graphics for finite element data. ! ! Changes: ! ! 04 April 2003 - Modified the format of the element and node data files. ! Element files no longer include a count of the number of elements and ! the number of nodes per element. Both node and element files now ! allow internal comment lines. ! ! 02 October 2000 - color graphics for 6 node elements were failing ! because SMAX and SMIN weren't passed in. ! ! 02 October 2000 - Fixed a bug and a flaw in the Q6 contour line routine. ! ! 19 September 2000 - Added option to write out element and node data. ! ! 12 June 2000 - Extensive FORTRAN 90 update. ! ! 15 April 2000 - FORTRAN 90 updates: logical operators, LEN_TRIM, ! and PARAMETER statements. ! ! 15 October 1999 - Increased loop count in BUZZ. Added IERROR ! flag in TRANSFORM. ! ! Here's an odd fact: The mapping ! ! I: 1 2 3 4 5 6 ! X: 38 40 40 39 39 40 ! Y: 8.5E+00 7.5 8.5 8.5 8.25 8.25 ! ! produces the transform: ! ! Y = - ETA**2 + 8.5E+00 ! ! which in turn produces a singular jacobian along the line ETA = 0, ! since dY/dETA = dY/dXSI = 0 there. ! ! 15 September 1999 - Added divergence contours and colors. ! ! 23 July 1999 - Tried to standardize the contour line calls by ! creating a new routine CONTOUR_LINE, also dropping the feature ! whereby the contour lines were sometimes varied in color. ! I still need to use S_CONTOUR and C_CONTOUR in the color contour ! routines instead of what I'm doing now. ! ! 15 July 1999 - Modifying contour level calculations so that ! the contour levels are computed all at once, and can be input ! individually by the user or from a file. ! ! 02 April 1999 - Rewrote EXAMPLE routine to make it depend ! on the number of elements in the X and Y directions, rather ! than the number of nodes. ! ! 01 April 1999 - Corrections to the stream line code and the ! 4 node quadrilateral contour line codes. ! ! Modified: ! ! 04 April 2003 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: maxbou = 800 integer, parameter :: maxcontour = 50 integer, parameter :: maxnpe = 6 integer, parameter :: maxnx = 151 integer, parameter :: maxny = 51 integer, parameter :: maxobj = 39 ! ! Parameters that depend on primitive parameters. ! integer, parameter :: maxelm = 2 * ( maxnx - 1 ) * ( maxny - 1 ) integer, parameter :: maxnp = ( 2 * maxnx - 1 ) * ( 2 * maxny - 1 ) character ( len = 10 ) arrow real bval integer c_contour(maxcontour) 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 = 80 ) element_file_name character ( len = 2 ) eqn(3,maxnp) real etaref(maxnpe) character ( len = 80 ) filgrf character ( len = 80 ) filinp character ( len = 80 ) filetec character ( len = 20 ) filtyp real grace real gval integer i integer icmax integer icmin integer icolor(maxobj) integer idata integer ierror integer ifile integer ios integer iplot character isay integer itable integer itemp integer iunit 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 line(maxobj) integer nbound integer ncontour integer nelem logical nflag(maxnp) logical nflag0(maxnp) logical nflag1(maxnp) integer node(maxnpe,maxelm) character ( len = 80 ) node_file_name logical none integer np integer npe integer numel(maxnp) character ( len = 40 ) object(maxobj) logical ovrlay real p(maxnp) real rho(maxnp) real rval real s(maxnp) real s2(maxnp) real scalee real scalen real scalev real s_contour(maxcontour) logical s_eqi logical show(maxobj) real smax real smin real srange real t(maxnp) real temp character ( len = 40 ) title character ( len = 40 ) title2 real u(maxnp) real v(maxnp) real x1max real x1min real x2max real x2min real x4max real x4min real xc(maxnp) real xmax real xmin real xsiref(maxnpe) real xsmax real xsmin real xtmax real xtmin real y1max real y1min 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 ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4' write ( *, '(a)' ) ' FORTRAN90 version' ! ! Greetings! ! call hello ( maxbou, maxcontour, maxelm, maxnp, maxnpe, maxnx, maxny, maxobj ) ! ! Set initial values. ! call init ( arrow, c_contour, cp, delx, dely, dev, dudxn, dudyn, & dvdxn, dvdyn, echo, eflag, eflagu, eqn, etaref, node_file_name, & filgrf, filinp, filtyp, grace, icmax, icmin, icolor, & idata, ifile, iplot, itable, iwrite, jcmax, & jcmin, labelx, labely, lbar, line, maxcontour, maxelm, maxnp, & maxnpe, maxobj, nbound, ncontour, nelem, nflag, nflag0, & nflag1, node, np, npe, object, ovrlay, p, rho, s_contour, & scalee, scalen, scalev, show, smax, smin, title, & title2, u, v, x1max, x1min, x2max, x2min, x4max, x4min, & xc, 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' ) ! ! DO FOREVER: Get a command from the user. ! do write ( *, '(a)' ) ' ' write ( *, '(a)' ) '? ("H" for help)' read ( *, '(a)', iostat = ios ) command if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if write ( 17, '(a)' ) trim ( command ) if ( echo ) then write ( *, '(a)' ) trim ( command ) end if if ( command == ' ') then cycle end if ! ! To simplify parsing, remove all blanks if possible. ! if ( .not. s_eqi ( command(1:5), 'TITLE' ) ) then call s_blank_delete ( command ) end if ! ! ARROWS = : Choose solid or line. ! 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Input request:' write ( *, '(a)' ) ' Choose HOLLOW, LINE, SOLID for arrows.' read ( *, '(a)' ) arrow write ( 17, '(a)' ) arrow if ( echo ) then write ( *, '(a)' ) arrow end if end if write ( *, '(a)' ) ' The arrow option is now ARROW = ' // arrow ! ! B: ! BO: switch the boundary display option. ! else if ( s_eqi ( command, 'B' ) .or. s_eqi ( command(1:2), 'BO' ) ) then show(1) = .not. show(1) if ( show(1) ) then write ( *, '(a)' ) 'The boundary will be shown.' else write ( *, '(a)' ) '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 ( *, '(a)' ) 'The background will be shown.' else write ( *, '(a)' ) '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 ( *, '(a)' ) 'The color bar will be shown.' else write ( *, '(a)' ) 'The color bar will NOT be shown.' end if ! ! BH: bottom half. ! else if ( s_eqi ( command, 'BH' ) ) then ysmax = ysmin + 0.5E+00 * ( 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.5E+00 * ( xsmax - xsmin ) ysmax = ysmin + 0.5E+00 * ( 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.5E+00 * ( 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.5E+00 * ( xsmax - xsmin ) ysmax = ysmin + 0.5E+00 * ( 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Number Color Name' write ( *, '(a)' ) ' ' do i = 1, maxobj write ( *, '(i2,2x,i3,2x,a)' ) i, icolor(i), trim ( object(i) ) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter an object number, and a color number.' read ( *, *, iostat = ios ) itemp, jtemp if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if write ( 17, * ) itemp, jtemp if ( echo ) then write ( *, * ) itemp, jtemp end if if ( 1 <= itemp .and. itemp <= maxobj ) then icolor(itemp) = jtemp else write ( *, '(a)' ) '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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Please use the DEV command first!' cycle end if if ( s_eqi ( command(1:3), 'CC=' ) ) then read ( command(4:), *, iostat = ios ) itable if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Built in color tables include:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) '1 low black to high white.' write ( *, '(a)' ) '2 low white to high black.' write ( *, '(a)' ) '3 low blue to high yellow.' write ( *, '(a)' ) '4 low red, high blue, with bands between.' write ( *, * ) '5 low red, yellow, green, blue, high white.' write ( *, * ) '6 low white, blue, green, yellow, high red' write ( *, * ) '7 low blue to high red.' write ( *, * ) '8 linear table between 2 user colors.' write ( *, * ) '9 linear table between N user colors.' write ( *, * ) ' ' write ( *, * ) 'Enter a color table index between 1 and 9,' write ( *, * ) 'or 0 to enter a color table from a file.' read ( *, *, iostat = ios ) itable if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if write ( 17, * ) itable if ( echo ) then write ( *, * ) itable end if end if call color_table_choose ( dev, echo, filgrf, icmax, icmin, ierror, & iplot, itable, ovrlay, x1max, x1min, y1max, y1min ) if ( itable == 1 .or. itable == 2 ) 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(1:5), 'COLOR' ) ) then write ( *, '(a)' ) 'Enter the color index between 0 and 255' read ( *, *, iostat = ios ) i if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if write ( 17, * ) i if ( echo ) then write ( *, * ) i end if write ( *, '(a)' ) 'Enter(R,G,B)' write ( *, '(a)' ) 'Note: (0,0,0) is black, (1,1,1) is white!' read ( *, *, iostat = ios ) rval, gval, bval if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if write ( 17, * ) rval, gval, bval if ( echo ) then write ( *, * ) rval, gval, bval end if call setclr ( i, rval, gval, bval ) ! ! 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 ( *, '(a)' ) 'CP colors will be plotted.' else write ( *, '(a)' ) '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 color_box if ( x1min == x1max ) then x1min = 0.0E+00 x1max = 1.0E+00 end if if ( y1min == y1max ) then y1min = 0.0E+00 y1max = 1.0E+00 end if ierror = 0 call setwcd ( x1min, y1min, x1max, y1max, ierror ) call buzz ( dev, x1min, x1max, y1min, y1max ) ! ! 'DEV = ' Choose the graphics device. ! else if ( s_eqi ( command(1:3), 'DEV' ) ) then if ( dev /= ' ' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Error!' write ( *, '(a)' ) ' You have already chosen device ' // trim ( dev ) write ( *, '(a)' ) ' You may not change your mind!' cycle end if if ( s_eqi ( command(1:4), 'DEV=' ) ) then dev = command(5:) else write ( *, * ) ' ' write ( *, * ) 'Enter the graphics device desired.' write ( *, * ) ' ' write ( *, * ) 'Options include:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CGMB output to a CGM binary file.' write ( *, '(a)' ) 'PS output to a PostScript file.' write ( *, '(a)' ) 'XWS output to an X window screen.' read ( *, '(a)', iostat = ios ) dev if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if 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 ( *, '(a)' ) 'Output will be to a CGM binary file "display.cgm".' else if ( s_eqi ( dev, 'PS' ) ) then write ( *, '(a)' ) 'Output will be to a PostScript file "display.ps".' else if ( s_eqi ( dev, 'XWS' ) ) then write ( *, '(a)' ) 'Output will be to an X window screen.' else write ( *, * ) 'Your device '//dev//' was not recognized!' dev = ' ' end if ! ! DIV: show divergence contours ! else if ( s_eqi ( command, 'DIV' ) ) then show(38) = .not. show(38) if ( show(38) ) then write ( *, '(a)' ) 'Divergence contours will be shown.' else write ( *, '(a)' ) 'Divergence contours will NOT be shown.' end if ! ! DIVC: show divergence color plots. ! else if ( s_eqi ( command, 'DIVC' ) ) then show(39) = .not. show(39) if ( show(39) ) then write ( *, '(a)' ) 'Divergence colors will be shown.' else write ( *, '(a)' ) 'Divergence colors will NOT be shown.' end if ! ! 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 ( *, '(a)' ) 'Elements will be shown.' else write ( *, '(a)' ) '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 ( *, '(a)' ) 'Element colors will be shown.' else write ( *, '(a)' ) '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 ( *, '(a)' ) 'Element numbers will be shown.' else write ( *, '(a)' ) 'Element numbers will NOT be shown.' end if ! ! EXAMple: Define example problem. ! else if ( s_eqi ( command(1:4), 'EXAM' ) ) then if ( ifile > 0 ) then close ( unit = 2 ) write ( *, * ) 'OpnFil is closing the old file ' & // trim ( node_file_name ) end if node_file_name = ' ' write ( *, * ) ' ' write ( *, * ) 'Enter number of nodes per element (3, 4, 6):' read ( *, * ) npe write ( 17, * ) npe if ( echo ) then write ( *, * ) npe end if call example ( eqn, maxelm, maxnp, maxnpe, nelem, node, np, npe, p, & rho, u, v, xc, yc ) nflag(1:np) = .true. nflag0(1:np) = .true. nflag1(1:np) = .true. eflag(1:nelem) = .true. eflagu(1:nelem) = .true. call rsize ( 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 bound_set ( eqn, jbound, maxbou, maxelm, maxnpe, nbound, & nelem, node, np, npe ) idata = 1 ! ! 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 ! ! FILE = : set the name of the graphics output file. ! else if ( s_eqi ( command(1:4), 'FILE' ) ) then if ( iplot > 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - 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)', iostat = ios ) filgrf if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if write ( 17, '(a)') filgrf if ( echo ) then write ( *, '(a)' ) filgrf end if end if write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - 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 (JEFF or TECPLOT).' read ( *, '(a)', iostat = ios ) filtyp if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if write ( 17, '(a)' ) trim ( filtyp ) if ( echo ) then write ( *, '(a)' ) trim ( filtyp ) end if end if write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Note:' write ( *, * ) ' The input file type is set to ' // trim ( filtyp ) ! ! FRAME: switch the frame display option. ! else if ( s_eqi ( command, 'FRAME' ) ) then show(3) = .not. show(3) if ( show(3) ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Note:' write ( *, * ) ' A frame will be shown.' else write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Note:' write ( *, * ) ' A frame will NOT be shown.' end if ! ! 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' ) .or. s_eqi ( command, 'GRAPH' ) ) then call graph ( arrow, c_contour, cp, delx, dely, dev, dudxn, & dudyn, dvdxn, dvdyn, echo, eflag, eflagu, etaref, filgrf, & icmax, icmin, icolor, iplot, itable, iwork1, & iwork2, jbound, jcmax, jcmin, lbar, line, maxbou, maxcontour, & maxelm, maxnp, maxnpe, maxobj, nbound, ncontour, nelem, & nflag, nflag0, nflag1, node, np, npe, numel, object, ovrlay, & p, rho, s, s_contour, s2, scalee, scalen, scalev, show, smax, & smin, srange, t, & title, title2, u, v, x1max, x1min, x2max, x2min, xc, 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:), *, iostat = ios ) grace if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if else write ( *, * ) 'Enter the grace margin:' read ( *, * ) grace write ( 17, *)grace if ( echo ) then write ( *, * ) grace end if end if write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - 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, maxcontour, maxelm, maxnp, maxnpe, maxnx, maxny, & maxobj ) ! ! 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:5), 'ICMAX' ) ) then if ( s_eqi ( command(1:6), 'ICMAX=' ) ) then read ( command(7:), *, iostat = ios ) icmax if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if else write ( *, '(a)' ) 'Enter ICMAX, maximum color index.' read ( *, * ) icmax end if if ( icmax > 255 ) then write ( *, '(a)' ) '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:5), 'ICMIN' ) ) then if ( s_eqi ( command(1:6), 'ICMIN=' ) ) then read ( command(7:), *, iostat = ios ) icmin if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if else write ( *, '(a)' ) 'Enter ICMIN, minimum color index.' read ( *, * ) icmin end if if ( icmin < 2 ) then write ( *, '(a)' ) 'ICMIN must be at least 2.' icmin = 2 end if write ( *, * ) 'Minimum available color set to ', icmin ! ! INIT: set variables to initial (zero!) values. ! else if ( s_eqi ( command(1:4), 'INIT' ) ) then call init ( arrow, c_contour, cp, delx, dely, dev, dudxn, dudyn, & dvdxn, dvdyn, echo, eflag, eflagu, eqn, etaref, node_file_name, & filgrf, filinp, filtyp, grace, icmax, icmin, icolor, & idata, ifile, iplot, itable, iwrite, jcmax, & jcmin, labelx, labely, lbar, line, maxcontour, maxelm, maxnp, & maxnpe, maxobj, nbound, ncontour, nelem, nflag, nflag0, & nflag1, node, np, npe, object, ovrlay, p, rho, s_contour, & scalee, scalen, scalev, show, smax, smin, title, & title2, u, v, x1max, x1min, x2max, x2min, x4max, x4min, & xc, xsiref, xsmax, xsmin, y1max, y1min, y2max, y2min, y4max, & y4min, yc, ysmax, ysmin ) ! ! IWRITE = : set the debugging output level. ! else if ( s_eqi ( command(1:6), 'IWRITE' ) ) then if ( s_eqi ( command(1:7), 'IWRITE=' ) ) then read ( command(8:), *, iostat = ios ) iwrite if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if else write ( *, '(a)' ) 'Enter IWRITE, debug output level.' read ( *, * ) iwrite end if 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:), *, iostat = ios ) jcmax if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Fatal error!' write ( *, '(a)' ) ' Error or end of file reading user input!' cycle end if 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:), *, iostat = ios ) jcmin if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if 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 velocity vectors will be shown.' else write ( *, * ) 'Kinematic velocity vectors 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 ! ! LIST: list current values ! else if ( s_eqi ( command(1:4), 'LIST' ) ) then call list ( delx, dely, dev, echo, node_file_name, grace, icmax, icmin, & icolor, idata, ifile, iplot, itable, iwrite, & jbound, maxbou, maxnp, maxobj, nbound, ncontour, nelem, np, & npe, object, p, scalev, show, title, title2, u, v, x2max, x2min, & xmax, xmin, xsmax, xsmin, ymax, & ymin, y2max, y2min, ysmax, ysmin ) ! ! LH: left half. ! else if ( s_eqi ( command, 'LH' ) ) then xsmax = xsmin + 0.5E+00 * ( xsmax - xsmin ) temp = 0.25E+00 * ( 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 ( *, '(i2,2x,i3,2x,a)' ) i, line(i), trim ( 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 ( *, *, iostat = ios ) itemp, jtemp if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if 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.25E+00 * ( xsmax - xsmin ) xsmin = xsmin + temp xsmax = xsmax - temp temp = 0.25E+00 * ( 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.25E+00 * ( 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.5E+00 * ( xsmax - xsmin ) temp = 0.25E+00 * ( 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.5E+00 * ( xsmax - xsmin ) temp = 0.25E+00 * ( 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' ) .or. s_eqi ( command, 'NODES' ) ) then show(4) = .not. show(4) if ( show(4) ) then write ( *, * ) 'Nodes will be shown.' else write ( *, * ) 'Nodes will NOT be shown.' end if ! ! NCONTOUR = : Set the number of contour lines. ! else if ( s_eqi ( command(1:2), 'NC' ) ) then if ( s_eqi ( command, 'NCONTOUR=' ) .or. s_eqi ( command, 'NCON=' ) ) then read ( command(6:), *, iostat = ios ) ncontour if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if write ( *, * ) 'Number of contour lines set to ',ncontour else write ( *, * ) 'Enter number of contour lines.' read ( *, * ) ncontour write ( 17, * ) ncontour if ( echo ) then write ( *, * ) ncontour 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 ! ! OVERLAY: Switch the overlay value. ! else if ( s_eqi ( command, 'OV' ) .or. 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' ) .or. s_eqi ( command, 'PRESSURE' ) ) 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 ! ! Q: quit. ! QUIT or QY: QUIT NOW! ! else if ( s_eqi ( command(1:1), 'Q' ) ) then if ( s_eqi ( command(2:2), 'Y' ) .or. s_eqi ( command, 'QUIT' ) ) 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(1:1), 'Y' ) ) then if ( iplot > 0 ) then call grfcls end if if ( ifile > 0 ) then close ( unit = 2 ) end if ! if ( .not. s_eqi ( dev, 'CGMB' ) ) then if ( .false. ) then call get_unit ( iunit ) open ( unit = iunit, file = 'cgmout', status = 'old', & iostat = ios ) if ( ios /= 0 ) then close ( unit = iunit, status = 'delete' ) end if end if close ( unit = 17 ) exit end if ! ! READELEMENT filename ! specify the element file to be read. ! else if ( s_eqi ( command(1:11), 'READELEMENT' ) ) then filtyp = 'jeff' if ( command(12:12) == '=' ) then element_file_name = command(13:) else if ( 0 < len_trim ( command(12:) ) ) then element_file_name = command(12:) else write ( *, * ) 'Enter the name of the element data file.' read ( *, '(a)' ) element_file_name write ( 17, '(a)' ) trim ( element_file_name ) if ( echo ) then write ( *, '(a)' ) trim ( element_file_name ) end if end if ! ! Read the number of elements. ! call file_line_count ( element_file_name, nelem ) if ( nelem <= 0 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Error!' write ( *, * ) ' This problem has zero elements!' cycle end if if ( nelem > maxelm ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Error!' write ( *, * ) ' This problem has too many elements!' write ( *, * ) ' Number of elements NELEM = ', nelem write ( *, * ) ' DISPLAY4 can handle up to MAXELM = ', maxelm cycle end if ! ! Read the number of nodes per element. ! call file_column_count ( element_file_name, npe ) if ( npe /= 3 .and. npe /= 4 .and. npe /= 6 ) then ierror = 1 write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Error!' write ( *, * ) ' Legal values of NPE are 3, 4 and 6.' write ( *, * ) ' User input value is NPE = ', npe cycle end if call read_element ( element_file_name, ierror, maxelm, maxnp, & maxnpe, nelem, node, np, npe ) if ( ierror /= 0 ) then ierror = 0 cycle end if ifile = 2 xsmin = 0.0E+00 xsmax = 0.0E+00 ysmin = 0.0E+00 ysmax = 0.0E+00 ! ! READNODE filename ! else if ( s_eqi ( command(1:8), 'READNODE' ) ) then if ( nelem <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4 - Error:' write ( *, '(a)' ) ' No elements have been defined yet.' write ( *, '(a)' ) ' Please enter the READ ELEMENT command first!' cycle end if if ( command(9:9) == '=' ) then node_file_name = command(10:) else if ( 0 < len_trim ( command(9:) ) ) then node_file_name = command(9:) else write ( *, * ) 'Enter the name of the node data file.' read ( *, '(a)' ) node_file_name write ( 17, '(a)' ) trim ( node_file_name ) if ( echo ) then write ( *, '(a)' ) trim ( node_file_name ) end if end if call read_node ( node_file_name, np, p, rho, u, v, xc, yc ) jfile = 2 call adjeff ( delx, dely, eflagu, eflag, eqn, grace, & ifile, jbound, maxbou, maxelm, maxnp, maxnpe, nbound, & nelem, nflag, nflag0, node, np, npe, srange, & x1max, x1min, x2max, x2min, xc, xmax, xmin, xsmax, & xsmin, xtmax, xtmin, y1max, y1min, y2max, y2min, yc, ymax, & ymin, ysmax, ysmin, ytmax, ytmin ) ! ! READTECPLOT filename ! specify the TECPLOT data file to be read. ! else if ( s_eqi ( command(1:11), 'READTECPLOT' ) ) then filtyp = 'tecplot' if ( command(12:12) == '=' ) then filetec = command(13:) else write ( *, * ) 'Enter the name of the TECPLOT data file.' read ( *, '(a)' ) filetec write ( 17, '(a)' ) filetec if ( echo ) then write ( *, '(a)' ) filetec end if end if call read_tecplot ( filetec, ierror, ifile, maxelm, maxnp, maxnpe, & nelem, node, np, npe, p, u, v, xc, yc ) if ( ierror /= 0 ) then ierror = 0 cycle end if xsmin = 0.0E+00 xsmax = 0.0E+00 ysmin = 0.0E+00 ysmax = 0.0E+00 ifile = 2 ! ! Initially, all elements will be visible. ! eflag(1:nelem) = .true. eflagu(1:nelem) = .true. nflag(1:np) = .true. nflag0(1:np) = .true. ! ! We give default values to the equations. ! eqn(1,1:np) = 'UM' eqn(2,1:np) = 'VM' eqn(3,1:np) = 'PC' ! ! Check element orientation. ! call element_check ( maxelm, maxnpe, nelem, node, np, npe, xc, yc ) ! ! Get the region size. ! call rsize ( 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 ) ! ! Determine the boundary edges. ! call bound_set ( eqn, jbound, maxbou, maxelm, maxnpe, nbound, & nelem, node, np, npe ) ! ! RH: right half. ! else if ( s_eqi ( command, 'RH') ) then temp = 0.50E+00 * ( 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, '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:), *, iostat = ios ) scalee if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if 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:), *, iostat = ios ) scalen if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if 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:), *, iostat = ios ) scalev if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if write ( *, * ) 'Velocity vector scale factor set to ', scalev ! ! TC: top center quarter. ! else if ( s_eqi ( command, 'TC' ) ) then temp = 0.25E+00 * ( xsmax - xsmin ) xsmin = xsmin + temp xsmax = xsmax - temp ysmin = ysmin + 0.5E+00 * ( 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.5E+00 * ( 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)', iostat = ios ) title2 if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if 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)', iostat = ios ) title if ( ios /= 0 ) then write ( *, * ) ' ' write ( *, * ) 'DISPLAY4 - Fatal error!' write ( *, * ) ' Error or end of file reading user input!' cycle end if 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.5E+00 * ( xsmax - xsmin ) ysmin = ysmin + 0.5E+00 * ( 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.5E+00 * ( xsmax - xsmin ) ysmin = ysmin + 0.5E+00 * ( ysmax - ysmin ) call pltbox ( grace, srange, x1max, x1min, x2max, x2min, xsmax, & xsmin, y1max, y1min, y2max, y2min, ysmax, ysmin ) ! ! 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 vizelm ( echo, 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 ( *, * ) 'DISPLAY4 - 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 contours. ! 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 ! ! WRITEELEMENT filename ! write element data to a file. ! else if ( s_eqi ( command(1:12), 'WRITEELEMENT' ) ) then if ( command(13:13) == '=' ) then element_file_name = command(14:) else if ( 0 < len_trim ( command(13:) ) ) then element_file_name = command(13:) else write ( *, * ) 'Enter the name of the element data file.' read ( *, '(a)' ) element_file_name write ( 17, '(a)' ) trim ( element_file_name ) if ( echo ) then write ( *, '(a)' ) trim ( element_file_name ) end if end if call write_element ( element_file_name, maxelm, maxnpe, nelem, node, npe ) ! ! WRITENODE filename ! write node data to a file. ! else if ( s_eqi ( command(1:9), 'WRITENODE' ) ) then if ( command(10:10) == '=' ) then node_file_name = command(11:) else if ( 0 < len_trim ( command(10:) ) ) then node_file_name = command(10:) else write ( *, * ) 'Enter the name of the output node data file.' read ( *, '(a)' ) node_file_name write ( 17, '(a)' ) trim ( node_file_name ) if ( echo ) then write ( *, '(a)' ) trim ( node_file_name ) end if end if call write_node ( node_file_name, np, p, u, v, xc, yc ) ! ! WRITE TECPLOT filename ! write data to a TECPLOT file. ! else if ( s_eqi ( command(1:12), 'WRITETECPLOT' ) ) then call write_tecplot ( maxelm, maxnpe, nelem, node, np, npe, p, rho, & u, v, xc, yc ) ! ! X: set the data window. ! else if ( s_eqi ( command, 'X' ) ) then call getwin ( 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 ! ! 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 ( *, * ) 'DISPLAY4 - Warning!' write ( *, * ) ' Your command was not recognized:' write ( *, '(a)' ) ' "' // trim ( command ) // '"' end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DISPLAY4:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine adjeff ( delx, dely, eflag, eflagu, eqn, grace, ifile, jbound, & maxbou, maxelm, maxnp, maxnpe, nbound, nelem, nflag, nflag0, node, np, & npe, srange, x1max, x1min, x2max, x2min, xc, xmax, xmin, xsmax, xsmin, & xtmax, xtmin, y1max, y1min, y2max, y2min,yc, ymax, ymin, ysmax, ysmin, & ytmax, ytmin ) !*****************************************************************************80 ! !! ADJEFF reads and checks the node data. ! ! Modified: ! ! 30 March 1999 ! ! 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), element visibility flags. ! ! 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: ! ! 'UM' The horizontal momentum equation. ! 'UW' The condition U = 0 applied at a node on a fixed wall. ! 'VM' The vertical momentum equation. ! 'VW' The condition V = 0 applied at a node on a fixed wall. ! 'PC' The continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! ! Input, real GRACE, the size of the "grace" margin. ! ! Input/output, integer IFILE. ! Records the status of the data file whose name is NODE_FILE_NAME. ! ! 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. ! ! Input, integer MAXBOU. ! The amount of storage available for the IBOUND array. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNP, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer MAXNY, the maximum allowed value of NY. ! ! Output, integer NBOUND, the number of points defining the boundary. ! ! Output, integer NELEM, the number of elements. ! ! Output, logical NFLAG(MAXNP), flags nodes which are active. ! ! Output, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Output, integer NP, the number of nodes. ! ! Output, integer NPE, the number of nodes per element. ! ! Output, real SRANGE. ! The maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! Output, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! 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. ! ! Output, real XC(MAXNP), the X coordinates of the nodes. ! ! Output, real XMAX. ! The maximum X coordinate of all the nodes. ! The maximum entry in the XC array. ! ! Output, real XMIN. ! The minimum X coordinate of all the nodes. ! The minimum entry in the XC array. ! ! 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. ! ! 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. ! ! Output, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! ! 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. ! ! Output, real YC(MAXNP), the Y coordinates of the nodes. ! ! Output, real YMAX, the maximum Y coordinate of all the nodes. ! ! Output, real YMIN, the minimum Y coordinate of all the nodes. ! ! 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. ! ! Output, real YSMIN, the minimum displayed Y coordinate. ! implicit none integer maxbou integer maxelm integer maxnp integer maxnpe real delx real dely logical eflag(maxelm) logical eflagu(maxelm) character ( len = 2 ) eqn(3,maxnp) real grace integer i integer ifile integer jbound(5,maxbou) integer nbound integer nelem logical nflag(maxnp) logical nflag0(maxnp) integer node(maxnpe,maxelm) integer np integer npe real srange 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 ifile = 2 ! ! Initially, all elements will be visible. ! eflag(1:nelem) = .true. eflagu(1:nelem) = .true. nflag(1:np) = .true. nflag0(1:np) = .true. ! ! We give default values to the equations. ! eqn(1,1:np) = 'UM' eqn(2,1:np) = 'VM' eqn(3,1:np) = 'PC' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ADJEFF - Note:' write ( *, '(a)' ) ' The node data has been read.' ! ! Check element orientation. ! call element_check ( maxelm, maxnpe, nelem, node, np, npe, xc, yc ) ! ! Get the region size. ! call rsize ( 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 ) ! ! Determine the boundary edges. ! call bound_set ( eqn, jbound, maxbou, maxelm, maxnpe, nbound, nelem, node, & np, npe ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ADJEFF - Read input data.' write ( *, * ) ' NP = ', np write ( *, * ) ' NPE = ', npe write ( *, * ) ' NELEM = ', nelem return end subroutine arrow_line ( xstart, ystart, xtip, ytip ) !*****************************************************************************80 ! !! ARROW_LINE makes a line drawing of an arrow at any point on a graph. ! ! Discussion: ! ! 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: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real XSTART, YSTART, the starting point for the arrow. ! ! Input, real XTIP, YTIP, the end point for the arrow. ! implicit none real alpha real del real dist real pi real theta real xbase real xleft real xrite real xstart real xtip real ybase real yleft real yrite real ystart real ytip ! ! left ! |\ ! | \ ! start ------- base tip ! | / ! |/ ! rite ! if ( xstart == xtip .and. ystart == ytip ) then return end if theta = 0.5E+00 * pi() - atan2 ( 2.0E+00, 1.0E+00 ) dist = sqrt ( ( xtip - xstart )**2 + ( ytip - ystart )**2 ) alpha = atan2 ( ytip - ystart, xtip - xstart ) del = sqrt ( 5.0E+00 ) * dist / 3.0E+00 call movcgm ( xstart, ystart ) xbase = xstart + dist * cos ( alpha ) * 2.0E+00 / 3.0E+00 ybase = ystart + dist * sin ( alpha ) * 2.0E+00 / 3.0E+00 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 arrow_poly ( arrow, xstart, ystart, xtip, ytip ) !*****************************************************************************80 ! !! ARROW_POLY makes a polygonal drawing of an arrow at any point on a graph. ! ! Discussion: ! ! 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: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) ARROW, specifies how the arrow will be drawn. ! LINE, line drawing; ! SOLID, polygonal shape, filled in; ! HOLLOW, polygonal shape, outlined. ! ! Input, real XSTART, YSTART, the starting point for the arrow. ! ! Input, real XTIP, YTIP, the end point for the arrow. ! implicit none integer, parameter :: npts = 7 real alpha character ( len = * ) arrow real dist real pi logical s_eqi real theta real xpts(npts+1) real xstart real xtip real ypts(npts+1) real ystart real ytip ! ! left ! 5 ! |\ ! 7--------------6 \ ! start base 4tip ! 1--------------2 / ! |/ ! 3 ! rite ! if ( xstart == xtip .and. ystart == ytip ) then return end if dist = sqrt ( ( xtip - xstart )**2 + ( ytip - ystart )**2 ) theta = 0.5E+00 * pi() - atan2 ( 2.0, 1.0E+00 ) alpha = atan2 ( ytip - ystart, xtip - xstart ) xpts(1) = xstart + dist * sin ( alpha ) / 10.0E+00 ypts(1) = ystart - dist * cos ( alpha ) / 10.0E+00 xpts(2) = xstart + dist * sin ( alpha ) / 10.0E+00 & + dist * cos ( alpha ) * 2.0E+00 / 3.0E+00 ypts(2) = ystart - dist * cos ( alpha ) / 10.0E+00 & + dist * sin ( alpha ) * 2.0E+00 / 3.0E+00 xpts(3) = xstart + dist * sqrt ( 5.0E+00 ) * cos ( alpha - theta ) / 3.0E+00 ypts(3) = ystart + dist * sqrt ( 5.0E+00 ) * sin ( alpha - theta ) / 3.0E+00 xpts(4) = xtip ypts(4) = ytip xpts(5) = xstart + dist * sqrt(5.0)*cos(alpha+theta)/3.0E+00 ypts(5) = ystart + dist * sqrt(5.0)*sin(alpha+theta)/3.0E+00 xpts(6) = xstart - dist * sin(alpha) / 10.0E+00 + dist * cos(alpha)*2.0/3.0E+00 ypts(6) = ystart + dist * cos(alpha) / 10.0E+00 + dist * sin(alpha)*2.0/3.0E+00 xpts(7) = xstart - dist * sin(alpha) / 10.0E+00 ypts(7) = ystart + dist * cos(alpha) / 10.0E+00 xpts(8) = xpts(1) ypts(8) = ypts(1) if ( s_eqi ( arrow, 'HOLLOW' ) ) then call plylin ( npts+1, xpts, ypts ) else if ( s_eqi ( arrow, 'SOLID' ) ) then call plygon ( npts, xpts, ypts ) end if return end subroutine bound_print ( jbound, maxbou, nbound ) !*****************************************************************************80 ! !! BOUND_PRINT prints the edges which are boundaries of the region. ! ! Modified: ! ! 05 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! 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. ! ! Input, integer MAXBOU, the amount of storage in the IBOUND array. ! ! Input, integer NBOUND, the number of points defining the boundary. ! implicit none integer maxbou integer i integer j integer jbound(5,maxbou) integer nbound write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BOUND_PRINT' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I, Element, N1, N2, N3, Wall/Open' write ( *, '(a)' ) ' ' do i = 1, nbound write ( *, * ) i, ( jbound(j,i), j = 1, 5 ) end do return end subroutine bound_set ( eqn, jbound, maxbou, maxelm, maxnpe, nbound, nelem, & node, np, npe ) !*****************************************************************************80 ! !! BOUND_SET finds edges which are boundaries of the region. ! ! Discussion: ! ! 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, BOUND_SET is bound to be slow. If you know ! where your boundaries occur, then you might want to use ! a simpler routine. However, BOUND_SET 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: ! ! 06 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! 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: ! ! 'UM' The horizontal momentum equation. ! 'UW' The condition U = 0 applied at a node on a fixed wall. ! 'VM' The vertical momentum equation. ! 'VW' The condition V = 0 applied at a node on a fixed wall. ! 'PC' The continuity equation. ! 'PB' The condition P = 0 applied at (XMAX,YMAX). ! ! 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. ! ! Input, integer MAXBOU. ! The amount of storage available for the IBOUND array. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NBOUND, the number of points defining the boundary. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NPE, the number of nodes per element. ! implicit none integer maxbou integer maxelm integer maxnpe integer np character ( len = 2 ) eqn(3,np) integer ielem 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 nelem integer node(maxnpe,maxelm) integer npe ! ! NEDGE is the number of sides an element has. ! 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 IELEM... ! do ielem = 1, nelem ! ! ...with edge J... ! do j = 1, nedge ! ! ...with corner nodes (N1, N2); ! n1 = node(j,ielem) if ( j < nedge ) then jp1 = j + 1 else jp1 = 1 end if n2 = node(jp1,ielem) ! ! compare with each element K not equal to I... ! do k = 1, nelem if ( k /= ielem ) then ! ! ...with edge L... ! do l = 1, nedge ! ! ...with corner nodes (M1,M2)... ! m1 = node(l,k) if ( l < nedge ) then lp1 = l+1 else lp1 = 1 end if 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BOUND_SET - Warning!' write ( *, '(a)' ) ' 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) = ielem jbound(2,nbound) = j if ( npe == 3 ) then jbound(3,nbound) = jp1 jbound(4,nbound) = 0 else if ( npe == 4 ) then jbound(3,nbound) = jp1 jbound(4,nbound) = 0 else if ( npe == 6 ) then if ( j == 1 ) then jbound(3,nbound) = 5 else if ( j == 2 ) then jbound(3,nbound) = 6 else if ( j == 3 ) then jbound(3,nbound) = 4 end if 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 ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BOUND_SET - Note:' write ( *, * ) ' Found NBOUND = ', nbound, ' boundary edges.' return end subroutine bound_show ( eflag, etaref, jbound, line, maxbou, maxnpe, maxobj, & nbound, nelem, nflag, node, np, npe, xc, xsiref, yc ) !*****************************************************************************80 ! !! BOUND_SHOW displays the boundary. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical EFLAG(NELEM), element visibility flags. ! ! Input, real ETAREF(MAXNPE), the ETA coordinates of the reference nodes. ! ! 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. ! ! Input, integer MAXBOU, the storage available for the IBOUND array. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NBOUND, the number of points defining the boundary. ! ! Input, integer NELEM, the number of elements. ! ! Input, logical NFLAG(NP), flags nodes which are active. ! ! Input, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NPE, the number of nodes per element. ! ! Input, real XC(NP), the X coordinates of the nodes. ! ! Input, real XSIREF(MAXNPE), the XSI coordinates of the reference nodes. ! ! Input, real YC(NP), the Y coordinates of the nodes. ! implicit none integer maxbou integer maxnpe integer maxobj integer nelem integer np logical eflag(nelem) real etaref(maxnpe) integer i integer ielem 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 "simple". ! if ( npe == 3 .or. npe == 4 ) 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 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 iso_line_q6 ( x1, y1, x2, y2, ielem, maxnpe, & nelem, node, np, npe, xc, yc ) else if ( jbound(5,i) == 0 ) then call iso_line_q6 ( 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 iso_line_q6 ( 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 iso_line_q6 ( 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 ) !*****************************************************************************80 ! !! 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. ! implicit none 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 ) !*****************************************************************************80 ! !! 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. ! implicit none 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 ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none character c integer itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end subroutine color_bar ( c_contour, icolor, maxcontour, maxobj, ncontour, & s_contour, srange, x1, x2, y1, y2 ) !*****************************************************************************80 ! !! COLOR_BAR draws a color bar in a given rectangle. ! ! Modified: ! ! 16 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ICOLOR(MAXOBJ), the color index for each object. ! ! Input, integer MAXOBJ, the number of graphical "objects". ! ! Input, integer NCONTOUR, the number of color contour ! regions drawn, and hence, the number of colors ! to be displayed in the color bar. ! ! 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. ! ! Input, real SRANGE, the maximum of XSMAX-XSMIN and YSMAX-YSMIN. ! This gives the size of a square containing the data ! window. ! ! Input, real X1, X2, Y1, Y2, specify the minimum and ! maximum X and Y coordinates of the color bar. ! implicit none integer maxcontour integer maxobj real angle integer c_contour(maxcontour) character chrrel*14 character chrtmp*14 real cwide character ( len = 6 ) flush integer i integer icolor(maxobj) integer ncontour real pwide real s_contour(maxcontour) 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 = 1, ncontour xl = ( real ( ncontour - i + 1 ) * x1 + real ( i - 1 ) * x2 ) / & real ( ncontour ) xr = ( real ( ncontour - i ) * x1 + real ( i ) * x2 ) / & real ( ncontour ) xcorn(1) = xl xcorn(2) = xr xcorn(3) = xr xcorn(4) = xl xcorn(5) = xl call filclr ( c_contour(i) ) 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.0E+00 chrtmp = chrrel ( s_contour(1) ) call s_blank_delete ( chrtmp ) angle = 0.0E+00 x = x1 y = y1 - 1.5E+00 * cwide pwide = srange flush = 'left' call s_plot ( angle, cwide, pwide, trim ( chrtmp ), x, y, flush ) chrtmp = chrrel ( s_contour(ncontour) ) call s_blank_delete ( chrtmp ) angle = 0.0E+00 x = x2 y = y1 - 1.5E+00 * cwide pwide = srange flush = 'right' call s_plot ( angle, cwide, pwide, trim ( chrtmp ), x, y, flush ) return end subroutine color_box !*****************************************************************************80 ! !! COLOR_BOX draws a 16 by 16 box of colors in the current color table. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: npoly = 5 real angle real cwide character ( len = 6 ) flush real grace integer i integer icolor integer ierror integer j real pwide character ( len = 3 ) string real x(npoly) real x1max real x1min real x2max real x2min real xtemp real y(npoly) real y1max real y1min real y2max real y2min real ytemp ! ! Set the coordinate system to be 0 < = X <= 16.0, ! and 0.0E+00 < = Y <= 16.0E+00 ! x2min = 0.0E+00 x2max = 16.0E+00 y2min = 0.0E+00 y2max = 16.0E+00 grace = 0.05 x1min = x2min - grace * ( x2max - x2min ) x1max = x2max + grace * ( x2max - x2min ) y1min = y2min - grace * ( y2max - y2min ) y1max = y2max + grace * ( y2max - y2min ) ierror = 0 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 ! ! Print numeric indices. ! icolor = 1 call linclr ( icolor ) icolor = 0 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) write ( string, '(i3)' ) icolor angle = 0.0E+00 cwide = 0.015 * ( x2max - x2min ) pwide = x2max - x2min xtemp = real ( j ) - 0.5E+00 ytemp = 16.5E+00 - real ( i ) flush = 'center' call s_plot ( angle, cwide, pwide, string, xtemp, ytemp, flush ) icolor = icolor + 1 call plylin ( 5, x, y ) end do end do return end subroutine chrcti2 ( string, intval, ierror, lchar ) !*****************************************************************************80 ! !! CHRCTI2 finds and reads an integer from a string. ! ! Discussion: ! ! CHRCTI2 is given a string which may contain one or more integers. ! Starting at the first character position, it looks for the first ! substring that could represent an integer. If it finds such a string, ! it returns the integer's value, and the position of the last character ! read. ! ! Examples: ! ! STRING INTVAL LCHAR ! ! 'Apollo 13' 13 9 ! ' 1 ' 1 6 ! '1A' 1 1 ! '12,34,56' 12 2 ! 'A1A2A3' 1 2 ! '-1E2ABCD' -1 2 ! '-X20ABCD' 20 4 ! '23.45' 23 2 ! ' N = 34, $' 34 7 ! 'Oops!' 0 0 ! ! Modified: ! ! 26 September 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) STRING, the string to be read. ! Reading will begin at position 1 and terminate at the end of the ! string, or when no more characters can be read to form a legal integer. ! Blanks, commas, or other nonnumeric data will, in particular, ! cause the conversion to halt. ! ! Output, integer INTVAL, the integer read from the string, or 0 ! if there was an error. ! ! Output, integer IERROR, 0 an integer was found, 1 if no integer found. ! ! Output, integer LCHAR, the last character of STRING that is part ! of the integer. ! implicit none character chrtmp integer i integer idig integer ierror integer ihave integer intval integer isgn integer iterm integer lchar integer nchar character ( len = * ) string nchar = len ( string ) ierror = 0 i = 0 isgn = 1 intval = 0 ihave = 0 iterm = 0 ! ! Examine the next character. ! 10 continue i = i + 1 if ( i > nchar ) then iterm = 1 else chrtmp = string(i:i) ! ! Minus sign. ! if ( chrtmp == '-' ) then if ( ihave == 0 ) then ihave = 1 isgn = -1 else iterm = 1 end if ! ! Plus sign. ! else if ( chrtmp == '+' ) then if ( ihave == 0 ) then ihave = 1 else iterm = 1 end if ! ! Digit. ! else if ( lge ( chrtmp, '0' ) .and. lle ( chrtmp, '9' ) ) then ihave = 2 call digten ( chrtmp, idig ) intval = 10 * intval + idig ! ! Blank or TAB. ! else if ( ihave == 2 ) then iterm = 1 else ihave = 0 end if end if end if if ( iterm /= 1 ) then go to 10 end if if ( ihave == 2 ) then lchar = i - 1 intval = isgn * intval else ierror = 0 lchar = 0 intval = 0 end if return end function chrrel ( rval ) !*****************************************************************************80 ! !! CHRREL represents a real using 14 right justified characters. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real RVAL, a real number. ! ! Output (through function value), character ( len = 14 ) CHRREL, ! a right-justified character variable containing the ! representation of RVAL, using a G14.6 format. ! implicit none 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 color_table_choose ( dev, echo, filgrf, icmax, icmin, ierror, & iplot, itable, ovrlay, x1max, x1min, y1max, y1min ) !*****************************************************************************80 ! !! COLOR_TABLE_CHOOSE gets the color table choice from the user. ! ! Modified: ! ! 21 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = 10 ) DEV, the graphics output device: ! cgmb - CGM binary file. ! ps - PostScript file. ! XWS - X window screen (interactive). ! ! Input, logical ECHO, is TRUE if user input should be echoed. ! ! Input, character ( len = 80 ) FILGRF, the name of the output ! graphics file. ! ! Input, integer ICMAX, ICMIN, the maximum and minimum color ! indices to use in color contour graphics. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IPLOT, the number of plots made so far. ! ! Output, integer ITABLE, the desired color table. ! ! 1: low black to high white ! 2: low blue to high yellow ! 3: low red, high blue, with bands between. ! 4: low red, yellow, green, blue, high white. ! 5: low white, blue, green, yellow, high red. ! 6: low blue to high red ! 7: linear table between 2 user colors. ! 8: linear table between N user colors. ! 9: low white to high black. ! ! Input, logical OVRLAY. ! If OVRLAY is true, then the next time that a plot is ! requested, a "new frame" command is suppressed, so that ! the new plot is shown on top of the previous one. ! ! Input, real X1MAX, X1MIN, the maximum and minimum X ! coordinates of the plot, which includes a small grace margin. ! ! Input, real Y1MAX, Y1MIN, the maximum and minimum Y ! coordinates of the plot, which includes a small grace margin. ! implicit none character ( len = 10 ) dev logical echo character ( len = 80 ) filcol character ( len = 80 ) filgrf integer icmax integer icmin integer ierror integer iplot integer itable logical s_eqi logical ovrlay real x1max real x1min real y1max real y1min if ( itable == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COLOR_TABLE_CHOOSE - Input request:' write ( *, '(a)' ) ' Enter name of color table file.' read ( *, '(a)', err = 90, end = 90 ) filcol write ( 17, '(a)' ) filcol if ( echo ) then write ( *, '(a)' ) filcol end if call getctb ( icmin, icmax, filcol, ierror ) end if call preplt ( dev, echo, filgrf, icmax, icmin, iplot, itable, ovrlay ) if ( iplot > 1 ) then call color_table_set ( echo, icmin, icmax, itable ) end if if ( s_eqi ( dev, 'XWS' ) ) then call color_box ! ! Maybe this will help? ! call preplt ( dev, echo, filgrf, icmax, icmin, iplot, itable, ovrlay ) call color_box ierror = 0 call setwcd ( x1min, y1min, x1max, y1max, ierror ) call buzz ( dev, x1min, x1max, y1min, y1max ) end if return 90 continue write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COLOR_TABLE_CHOOSE - Error' ierror = 1 return end subroutine color_table_set ( echo, icmin, icmax, itable ) !*****************************************************************************80 ! !! COLOR_TABLE_SET sets up the color table. ! ! Discussion: ! ! This routine replaces the unreliable DRAWCGM routine SETCTB. ! ! The routine sets the colors between ICMIN and ICMAX, which ! should typically be 2 and 255. ! ! It will also set the values of color 0 to white, and ! color 1 to black. ! ! Modified: ! ! 21 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical ECHO, is TRUE if user input should be echoed. ! ! Input, integer ICMIN, ICMAX, the minimum and maximum color indices. ! ! Input, integer ITABLE, the desired table. ! ! 1: low black to high white ! 2: low white to high black. ! 3: low blue to high yellow ! 4: low red, high blue, with bands between. ! 5: low red, yellow, green, blue, high white. ! 5: low white, blue, green, yellow, high red. ! 7: low blue to high red ! 8: linear table between 2 user colors. ! 9: linear table between n user colors. ! implicit none real bhi real blo real bval logical echo real ghi real glo real gval integer i integer icmax integer icmin integer icol1 integer icol2 integer itable integer ival real pi real rhi real rlo real rval real temp real theta icmin = max ( icmin, 2 ) icmax = min ( icmax, 255 ) ! ! 1: Low black to high white ! 2: Low white to high black ! 3: Low blue to high yellow. ! 4: Low red, high blue, with bands between. ! 5: Low red, yellow, green, blue, high white. ! 6: Low white, blue, green, yellow, high red. ! 7: Low blue to high red ! if ( itable >= 1 .and. itable <= 7 ) then do i = icmin, icmax if ( icmin == icmax ) then temp = 0.5E+00 else temp = real ( i - icmin ) / real ( icmax - icmin ) end if if ( itable == 1 ) then bval = temp gval = temp rval = temp else if ( itable == 2 ) then bval = 1.0E+00 - temp gval = 1.0E+00 - temp rval = 1.0E+00 - temp else if ( itable == 3 ) then rval = temp gval = temp bval = 1.0E+00 - temp else if ( itable == 4 ) then theta = 0.5E+00 * pi() * temp rval = cos ( theta )**2 bval = sin ( theta )**2 gval = 0.8 * sin ( 10.0E+00 * theta )**6 else if ( itable == 5 ) then theta = 4.0E+00 * temp rval = exp(-(theta-1.0)**2) + exp(-(theta-4.0)**2) gval = exp(-(theta-2.0)**2) + exp(-(theta-4.0)**2) bval = exp(-(theta-3.0)**2) + exp(-(theta-4.0)**2) rval = max ( rval, 0.0E+00 ) rval = min ( rval, 1.0E+00 ) gval = max ( gval, 0.0E+00 ) gval = min ( gval, 1.0E+00 ) bval = max ( bval, 0.0E+00 ) bval = min ( bval, 1.0E+00 ) else if ( itable == 6 ) then theta = 4.0E+00 * temp rval = exp(-(theta-1.0)**2) + exp(-(theta-4.0)**2) gval = exp(-(theta-2.0)**2) + exp(-(theta-4.0)**2) bval = exp(-(theta-3.0)**2) + exp(-(theta-4.0)**2) rval = max ( rval, 0.0E+00 ) rval = min ( rval, 1.0E+00 ) gval = max ( gval, 0.0E+00 ) gval = min ( gval, 1.0E+00 ) bval = max ( bval, 0.0E+00 ) bval = min ( bval, 1.0E+00 ) else if ( itable == 7 ) then rval = temp gval = 0.0E+00 bval = 1.0E+00 - temp end if ival = i call setclr ( ival, rval, gval, bval ) end do ! ! 8: Interpolate between two values. ! else if ( itable == 8 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COLOR_TABLE_SET - Input request:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Enter Rlo, Glo, Blo, Rhi, Ghi, Bhi' write ( *, '(a)' ) ' Note: 0,0,0 is black, and 1,1,1 is white!' write ( *, '(a)' ) ' ' read ( *, *, end = 1952, err = 1964 ) rlo, glo, blo, rhi, ghi, bhi write ( 17, * ) rlo, glo, blo, rhi, ghi, bhi if ( echo ) then write ( *, * ) rlo, glo, blo, rhi, ghi, bhi end if rlo = max ( rlo, 0.0E+00 ) rhi = min ( rhi, 1.0E+00 ) glo = max ( glo, 0.0E+00 ) ghi = min ( ghi, 1.0E+00 ) blo = max ( blo, 0.0E+00 ) bhi = min ( bhi, 1.0E+00 ) do i = icmin, icmax if ( icmin == icmax ) then temp = 0.5E+00 else temp = real ( i - icmin ) / real ( icmax - icmin ) end if rval = rlo * ( 1.0E+00 - temp ) + rhi * temp gval = glo * ( 1.0E+00 - temp ) + ghi * temp bval = blo * ( 1.0E+00 - temp ) + bhi * temp ival = i call setclr ( ival, rval, gval, bval ) end do ! ! 9: Interpolate between several values. ! else if ( itable == 9 ) then icol1 = icmin write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COLOR_TABLE_SET - Input request:' write ( *, '(a)' ) ' ' write ( *, * ) ' Enter (R, G, B) for color index ', icol1 write ( *, * ) ' (0, 0, 0) is black.' write ( *, * ) ' (1, 1, 1) is white.' read ( *, * ) rlo, glo, blo write ( 17, * ) rlo, glo, blo if ( echo ) then write ( *, * ) rlo, glo, blo end if rlo = max ( rlo, 0.0E+00 ) glo = max ( glo, 0.0E+00 ) blo = max ( blo, 0.0E+00 ) do do write ( *, * ) ' Enter index of next color to set' write ( *, * ) ' between ', icol1+1, ' and ', icmax read ( *, * ) icol2 write ( 17, * ) icol2 if ( echo ) then write ( *, * ) icol2 end if if ( icol2 > icol1 .and. icol2 <= icmax ) then exit end if write ( *, * ) ' ' write ( *, * ) 'SETTAB - Warning!' write ( *, * ) ' Your color index was not accepted!' end do write ( *, * ) ' ' write ( *, * ) 'Enter (R, G, B) for color index ', icol2 read ( *, * ) rhi, ghi, bhi write ( 17, * ) rhi, ghi, bhi if ( echo ) then write ( *, * ) rhi, ghi, bhi end if rhi = min ( rhi, 1.0E+00 ) ghi = min ( ghi, 1.0E+00 ) bhi = min ( bhi, 1.0E+00 ) do i = icol1, icol2 if ( icol1 == icol2 ) then temp = 0.5E+00 else temp = real ( i - icol1 ) / real ( icol2 - icol1 ) end if rval = rlo * ( 1.0E+00 - temp ) + rhi * temp gval = glo * ( 1.0E+00 - temp ) + ghi * temp bval = blo * ( 1.0E+00 - temp ) + bhi * temp ival = i call setclr ( ival, rval, gval, bval ) end do if ( icol2 >= icmax ) then exit end if icol1 = icol2 rlo = rhi glo = ghi blo = bhi end do ! ! Unknown table. ! else write ( *, * ) ' ' write ( *, * ) 'COLOR_TABLE_SET - Fatal error!' write ( *, * ) ' Legal color table indices are ' write ( *, * ) ' between 1 and 9. Your value was ', itable end if ! ! Background color 0 is to be white. ! ival = 0 rval = 1.0E+00 gval = 1.0E+00 bval = 1.0E+00 call setclr ( ival, rval, gval, bval ) ! ! Foreground color 1 is to be black. ! ival = 1 rval = 0.0E+00 gval = 0.0E+00 bval = 0.0E+00 call setclr ( ival, rval, gval, bval ) return 1952 continue write ( *, * ) ' ' write ( *, * ) 'COLOR_TABLE_SET - Fatal error!' write ( *, * ) ' Unexpected end of file!' stop 1964 continue write ( *, * ) ' ' write ( *, * ) 'COLOR_TABLE_SET - Warning!' write ( *, * ) ' Illegal format for input data!' return end subroutine contour_color_l3 ( c_contour, eflag, maxcontour, maxelm, maxnpe, & ncontour, nelem, node, np, s, s_contour, xc, yc ) !*****************************************************************************80 ! !! CONTOUR_COLOR_L3 does color contours in a 3 node linear element. ! ! Modified: ! ! 16 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical EFLAG(NELEM), element visibility flags. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NCONTOUR, the number of color contours to use. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Input, integer NP, the number of nodes. ! ! Input, real S(NP), a scalar quantity associated with the nodes. ! This routine only looks at the values associated with ! corner element nodes. ! ! Input, real XC(NP), the X coordinates of the nodes. ! ! Input, real YC(NP), the Y coordinates of the nodes. ! implicit none integer maxcontour integer maxelm integer maxnpe integer np integer c_contour(maxcontour) logical eflag(maxelm) integer i integer i1 integer i2 integer i3 integer ncontour integer nelem integer node(maxnpe,maxelm) real s(np) real s_contour(maxcontour) real xc(np) real yc(np) ! ! Consider element I: ! do i = 1, nelem ! ! If it is visible... ! if ( eflag(i) ) then i1 = node(1,i) i2 = node(2,i) i3 = node(3,i) call tric_l3 ( c_contour, i1, i2, i3, maxcontour, ncontour, np, s, & s_contour, xc, yc ) end if end do return end subroutine contour_color_l4 ( c_contour, eflag, maxcontour, maxelm, maxnpe, & ncontour, nelem, node, np, s, s_contour, xc, yc ) !*****************************************************************************80 ! !! CONTOUR_COLOR_L4 does color contours for linear quadrilateral elements. ! ! Modified: ! ! 16 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical EFLAG(NELEM), element visibility flags. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NCONTOUR, the number of color contours to use. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Input, integer NP, the number of nodes. ! ! Input, real S(NP), a scalar quantity associated with the nodes. ! This routine only looks at the values associated with ! corner element nodes. ! ! Input, real XC(NP), the X coordinates of the nodes. ! ! Input, real YC(NP), the Y coordinates of the nodes. ! implicit none integer maxcontour integer maxnpe integer maxelm integer np integer c_contour(maxcontour) logical eflag(maxelm) integer i integer i1 integer i2 integer i3 integer ncontour integer nelem integer node(maxnpe,maxelm) real s(np) real s_contour(maxcontour) 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 tric_l3 ( c_contour, i1, i2, i3, maxcontour, ncontour, np, s, & s_contour, xc, yc ) i1 = node(3,i) i2 = node(4,i) i3 = node(1,i) call tric_l3 ( c_contour, i1, i2, i3, maxcontour, ncontour, np, s, & s_contour, xc, yc ) end if end do return end subroutine contour_color_q6 ( c_contour, eflag, etaref, jcmax, jcmin, & maxcontour, maxelm, maxnpe, ncontour, nelem, node, np, npe, s, s_contour, & smax, smin, xc, xsiref, yc ) !*****************************************************************************80 ! !! CONTOUR_COLOR_Q6 uses color to indicate values greater than a given value. ! ! Discussion: ! ! This routine is used for quantities associated with the six ! corner nodes of a 6 node finite element, in particular, ! vorticity or velocity magnitude. ! ! The triangular element is broken up into four three ! node triangles for further treatment. ! ! 2 ! |\ ! | \ ! 6--5 ! |\ |\ ! | \| \ ! 3--4--1 ! ! Thus, the four triangles are defined by the nodes as follows: ! ! Triangle Nodes ! ! 1 1, 5, 4 ! 2 2, 6, 5 ! 3 3, 4, 6 ! 4 4, 5, 6 ! ! Modified: ! ! 06 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical EFLAG(NELEM), element visibility flags. ! ! Input, real ETAREF(MAXNPE), the ETA coordinates of the reference nodes. ! ! Input, integer JCMAX, JCMIN, the maximum and minimum color ! indices to use for contours. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NCONTOUR, the number of contours to draw. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Input, integer NP, the number of nodes. ! ! Input, real S(NP), the value of S at each node. ! ! Input, real SMAX, SMIN, the maximum and minimum values of S. ! ! Input, real XC(NP), the X coordinates of the nodes. ! ! Input, real XSIREF(MAXNPE), the XSI coordinates of the reference nodes. ! ! Input, real YC(NP), the Y coordinates of the nodes. ! implicit none integer maxcontour integer maxnpe integer maxelm integer np integer c_contour(maxcontour) logical eflag(maxelm) real eta1 real eta2 real eta3 real etaref(maxnpe) integer i integer i1 integer i2 integer i3 integer ielem integer j integer jcmax integer jcmin integer ncontour integer nelem integer node(maxnpe,maxelm) integer npe integer nsub(4,3) real s(np) real s_contour(maxcontour) 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) = 5 nsub(1,3) = 4 nsub(2,1) = 2 nsub(2,2) = 6 nsub(2,3) = 5 nsub(3,1) = 3 nsub(3,2) = 4 nsub(3,3) = 6 nsub(4,1) = 4 nsub(4,2) = 5 nsub(4,3) = 6 ! ! Consider each element. ! do i = 1, nelem ielem = i ! ! ...proceed if the element is visible. ! if ( eflag(i) ) then ! ! Consider each of the four subtriangles. ! 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 tric_q6 ( c_contour, eta1, eta2, eta3, ielem, jcmax, jcmin, & maxcontour, maxnpe, ncontour, nelem, node, np, npe, & s_contour, s1, s2, s3, xc, xsi1, xsi2, xsi3, yc ) end do end if end do return end subroutine contour_levels ( c_contour, echo, jcmax, jcmin, maxcontour, name, & ncontour, s_contour, stmax, stmin, svmax, svmin ) !*****************************************************************************80 ! !! CONTOUR_LEVELS allows the user to set the contour levels. ! ! Modified: ! ! 16 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical ECHO, is TRUE if user input should be echoed. ! ! Input, character ( len = * ) NAME, the name of the variable. ! implicit none integer maxcontour integer c_contour(maxcontour) character ( len = 80 ) contour_filename integer contour_unit logical echo integer i character isay integer jcmax integer jcmin logical s_eqi character ( len = * ) name integer ncontour real s_contour(maxcontour) real stmax real stmin real sumax real sumin real svmax real svmin real x 10 continue write ( *, * ) ' ' write ( *, * ) name write ( *, * ) ' ' write ( *, * ) 'CONTOUR_LEVELS - Input request:' write ( *, * ) ' Total data range is ', stmin, ' to ', stmax write ( *, * ) ' Visible data range is ', svmin, ' to ', svmax write ( *, * ) ' ' write ( *, * ) ' F=File, T=total, U=user, V=visible' read ( *, '(a)', err = 40, end = 40 ) isay write ( 17, '(a)' ) isay if ( echo ) then write ( *, '(a)' ) isay end if ! ! File: ! if ( s_eqi ( isay, 'F' ) ) then write ( *, * ) ' ' write ( *, * ) 'Enter the filename containing contour values.' read ( *, '(a)' ) contour_filename call get_unit ( contour_unit ) open ( unit = contour_unit, file = contour_filename, status = 'old', & err = 50 ) ncontour = 0 20 continue if ( ncontour < maxcontour ) then read ( contour_unit, *, end = 30 ) x ncontour = ncontour + 1 s_contour(ncontour) = x go to 20 end if 30 continue close ( unit = contour_unit ) ! ! Total range limits: ! else if ( s_eqi ( isay, 'T' ) ) then do i = 1, ncontour s_contour(i) = ( real ( ncontour + 1 - i ) * stmin & + real ( i ) * stmax ) / real ( ncontour + 1 ) end do ! ! User input: ! else if ( s_eqi ( isay, 'U' ) ) then write ( *, * ) ' ' write ( *, * ) ' Enter new minimum and maximum values for contours.' read ( *, *, end = 40, err = 40 ) sumin, sumax write ( 17, * ) sumin, sumax if ( echo ) then write ( *, * ) sumin, sumax end if do i = 1, ncontour s_contour(i) = ( real ( ncontour + 1 - i ) * sumin & + real ( i ) * sumax ) / real ( ncontour + 1 ) end do ! ! Visible limits: ! else if ( s_eqi ( isay, 'V' ) ) then do i = 1, ncontour s_contour(i) = ( real ( ncontour + 1 - i ) * svmin & + real ( i ) * svmax ) / real ( ncontour + 1 ) end do ! ! Huh? ! else write ( *, * ) ' ' write ( *, * ) 'Your choice was not acceptable.' go to 10 end if ! ! Set the colors. ! do i = 1, ncontour c_contour(i) = ( ( ncontour + 1 - i ) * jcmin + i * jcmax ) & / ( ncontour + 1 ) end do return 40 continue write ( *, * ) ' ' write ( *, * ) 'CONTOUR_LEVELS - Warning!' write ( *, * ) ' There was trouble reading your input.' go to 10 50 continue write ( *, * ) ' ' write ( *, * ) 'CONTOUR_LEVELS - Warning!' write ( *, * ) ' Could not open the file.' go to 10 end subroutine contour_line ( c_contour, echo, eflag, etaref, jcmax, jcmin, & jcolor, maxcontour, maxelm, maxnp, maxnpe, ncontour, nelem, nflag, node, & np, npe, object_name, s, s_contour, xc, xsiref, yc ) !*****************************************************************************80 ! !! CONTOUR_LINE supervises the drawing of a contour line plot. ! ! ! Modified: ! ! 13 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none integer maxcontour integer maxelm integer maxnp integer maxnpe integer c_contour(maxcontour) logical echo logical eflag(maxelm) real etaref(maxnpe) integer i integer jcmax integer jcmin integer jcolor integer ncontour integer nelem logical nflag(maxnp) integer node(maxnpe,maxelm) integer np integer npe character ( len = 40 ) object_name real s(maxnp) real s_contour(maxcontour) real scon real stmax real stmin real svmax real svmin real xc(maxnp) real xsiref(maxnpe) real yc(maxnp) call linclr ( jcolor ) call fsize ( nflag, np, s, stmax, stmin, svmax, svmin ) if ( stmax <= stmin ) then write ( *, * ) ' ' write ( *, * ) 'GRAPH - Warning!' write ( *, * ) ' Could not display ' // trim ( object_name ) write ( *, * ) ' All data values were equal to ', stmin return end if call contour_levels ( c_contour, echo, jcmax, jcmin, maxcontour, & object_name, ncontour, s_contour, stmax, stmin, svmax, svmin ) do i = 1, ncontour scon = s_contour(i) if ( npe == 3 ) then call contour_line_l3 ( eflag, maxelm, maxnpe, nelem, node, & np, s, scon, xc, yc ) else if ( npe == 4 ) then call contour_line_l4 ( eflag, maxelm, maxnpe, nelem, node, & np, s, scon, xc, yc ) else if ( npe == 6 ) then call contour_line_q6 ( eflag, etaref, maxnpe, nelem, node, & np, npe, s, scon, xc, xsiref, yc ) end if end do return end subroutine contour_line_l3 ( eflag, maxelm, maxnpe, nelem, node, np, s, scon, & xc, yc ) !*****************************************************************************80 ! !! CONTOUR_LINE_L3 draws a single contour line in a single 3 node element. ! ! Modified: ! ! 06 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical EFLAG(NELEM), element visibility flags. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Input, integer NP, the number of nodes. ! ! Input, real S(NP), the quantity for which a contour line is desired. ! ! Input, real SCON, the value associated with the contour line. ! ! Input, real XC(MAXNP), the X coordinates of nodes. ! ! Input, real YC(MAXNP), the Y coordinates of nodes. ! implicit none integer maxnpe integer maxelm integer np logical eflag(maxelm) integer ielem integer n_cross integer n1 integer n2 integer n3 integer nelem integer node(maxnpe,maxelm) real s(np) real s1 real s2 real s3 real scon real smax real smin real x1 real x2 real x3 real xc(np) real y1 real y2 real y3 real yc(np) n_cross = 0 ! ! Draw the contour line by searching over each element. ! do ielem = 1, nelem if ( eflag(ielem) ) then n1 = node(1,ielem) x1 = xc(n1) y1 = yc(n1) s1 = s(n1) n2 = node(2,ielem) x2 = xc(n2) y2 = yc(n2) s2 = s(n2) n3 = node(3,ielem) x3 = xc(n3) y3 = yc(n3) s3 = s(n3) smin = min ( s1, s2, s3 ) smax = max ( s1, s2, s3 ) if ( smin <= scon .and. scon <= smax ) then n_cross = n_cross + 1 call cross_l3 ( s1, s2, s3, scon, x1, x2, x3, y1, y2, y3 ) end if end if end do write ( *, * ) 'N_CROSS = ', n_cross return end subroutine contour_line_l4 ( eflag, maxelm, maxnpe, nelem, node, np, s, scon, & xc, yc ) !*****************************************************************************80 ! !! CONTOUR_LINE_L4 draws a contour line in a 4 node rectangular element. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical EFLAG(NELEM), element visibility flags. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Input, integer NP, the number of nodes. ! ! Input, real S(NP), a scalar quantity associated with the nodes. ! ! Input, real SCON, the value associated with the contour line. ! ! Input, real XC(NP), the X coordinates of the nodes. ! ! Input, real YC(NP), the Y coordinates of the nodes. ! implicit none integer maxnpe integer maxelm integer np logical eflag(maxelm) integer i integer n1 integer n2 integer n3 integer n4 integer nelem integer node(maxnpe,nelem) real s(np) real s1 real s2 real s3 real s4 real scon real smax real smin real x1 real x2 real x3 real x4 real xc(np) real y1 real y2 real y3 real y4 real yc(np) ! ! Look at a particular element, IELEM. ! do i = 1, nelem ! ! Is the element visible? ! if ( eflag(i) ) then n1 = node(1,i) x1 = xc(n1) y1 = yc(n1) s1 = s(n1) smin = s1 smax = s1 n2 = node(2,i) x2 = xc(n2) y2 = yc(n2) s2 = s(n2) smin = min ( smin, s2 ) smax = max ( smax, s2 ) n3 = node(3,i) x3 = xc(n3) y3 = yc(n3) s3 = s(n3) smin = min ( smin, s3 ) smax = max ( smax, s3 ) n4 = node(4,i) x4 = xc(n4) y4 = yc(n4) s4 = s(n4) smin = min ( smin, s4 ) smax = max ( smax, s4 ) if ( smin <= scon .and. scon <= smax ) then call cross_q4 ( s1, s2, s3, s4, scon, x1, x2, x3, x4, y1, y2, y3, y4 ) end if end if end do return end subroutine contour_line_q6 ( eflag, etaref, maxnpe, nelem, node, np, npe, s, & scon, xc, xsiref, yc ) !*****************************************************************************80 ! !! CONTOUR_LINE_Q6 draws a contour line in a 6 node quadratic triangular element. ! ! Note: ! ! THIS ROUTINE SHOULD BE REWRITTEN SO THAT THE NUMBER OF SUBTRIANGLES ! MAY BE EASILY VARIED...I JUST NEED TO DEFINE THE X, Y, S ! values of each subtriangle. ! ! Discussion: ! ! The original six node quadratic element is broken up into four ! three node elements, which are treated as though they were ! linear: ! ! 2 ! |\ ! | \ ! 6--5 ! |\ |\ ! | \| \ ! 3--4--1 ! ! Thus, the four triangles are defined by the nodes as follows: ! ! Triangle Nodes ! ! 1 1, 5, 4 ! 2 5, 2, 6 ! 3 4, 6, 3 ! 4 6, 4, 5 ! ! Modified: ! ! 02 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, logical EFLAG(NELEM), element visibility flags. ! ! Input, real ETAREF(MAXNPE), the ETA coordinates of the reference nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NELEM, the number of elements. ! ! Input, logical NFLAG(NP), is TRUE for nodes which are visible. ! ! Input, integer NODE(MAXNPE,MAXELM), the global node numbers in ! each element. ! ! Input, integer NP, the number of nodes. ! ! Input, integer NPE, the number of nodes per element. ! ! Input, real S(NP), a scalar quantity associated with the nodes. ! ! Input, real SCON, the value associated with the contour line. ! ! Input, real XC(NP), the X coordinates of the nodes. ! ! Input, real XSIREF(MAXNPE), the XSI coordinates of the reference nodes. ! ! Input, real YC(NP), the Y coordinates of the nodes. ! implicit none integer maxnpe integer nelem integer np logical eflag(nelem) real etaref(maxnpe) integer i integer icross integer ielem logical inside integer j integer k integer l integer lines(3,4,2) integer m1 integer m2 integer n1 integer n2 integer node(maxnpe,nelem) integer npe real s(np) real s1 real s2 real scon real x1 real x2 real xc(np) real xsiref(maxnpe) real xx(6) real xx2(2) real y1 real y2 real yc(np) real yy(6) real yy2(2) ! ! List the local node numbers, which define the three lines, ! which make up each of the four triangles which subdivide ! the original 6 node triangular element 3 node triangles. ! lines(1,1,1) = 1 lines(1,1,2) = 5 lines(2,1,1) = 5 lines(2,1,2) = 4 lines(3,1,1) = 4 lines(3,1,2) = 1 lines(1,2,1) = 5 lines(1,2,2) = 2 lines(2,2,1) = 2 lines(2,2,2) = 6 lines(3,2,1) = 6 lines(3,2,2) = 5 lines(1,3,1) = 4 lines(1,3,2) = 6 lines(2,3,1) = 6 lines(2,3,2) = 3 lines(3,3,1) = 3 lines(3,3,2) = 4 lines(1,4,1) = 6 lines(1,4,2) = 4 lines(2,4,1) = 4 lines(2,4,2) = 5 lines(3,4,1) = 5 lines(3,4,2) = 6 ! ! To draw the contour line: ! ! Consider each element. ! do i = 1, nelem ! ! ...but proceed only if the element is visible. ! if ( eflag(i) ) then ! ! Consider each of the four subtriangles that compose the triangle. ! do k = 1, 4 icross = 0 ! ! Consider each side of a subtriangle. ! do j = 1, 3 m1 = lines(j,k,1) n1 = node(m1,i) m2 = lines(j,k,2) n2 = node(m2,i) s1 = s(n1) s2 = s(n2) ! ! Does the value SCON lie between the values at the endpoints of ! a side of a subtriangle? ! call line_seg_contains_point_1d ( s1, scon, s2, inside ) if ( inside ) then x1 = xsiref(m1) y1 = etaref(m1) x2 = xsiref(m2) y2 = etaref(m2) ! ! Consider 4 cases: ! ! A) S1 = S2 = SCON, the whole side is an intersection. ! B) S1 = SCON, ! C) S2 = SCON, ! D) SCON strictly between S1 and S2. Find the location of the crossing. ! if ( s1 == scon .and. s2 == scon ) then icross = icross + 1 xx(icross) = x1 yy(icross) = y1 icross = icross + 1 xx(icross) = x2 yy(icross) = y2 else if ( s1 == scon ) then icross = icross + 1 xx(icross) = x1 yy(icross) = y1 else if ( s2 == scon ) then icross = icross + 1 xx(icross) = x2 yy(icross) = y2 else icross = icross + 1 xx(icross) = x1 + ( scon - s1 ) * ( x2 - x1 ) / ( s2 - s1 ) yy(icross) = y1 + ( scon - s1 ) * ( y2 - y1 ) / ( s2 - s1 ) end if end if end do ! ! How many times did the contour value cross the subtriangle? ! if ( icross > 1 ) then do j = 1, icross do l = j+1, icross ielem = i if ( npe == 3 ) then xx2(1) = xx(j) yy2(1) = yy(j) xx2(2) = xx(l) yy2(2) = yy(l) call plylin ( 2, xx2, yy2 ) else x1 = xx(j) y1 = yy(j) x2 = xx(l) y2 = yy(l) call iso_line_q6 ( x1, y1, x2, y2, ielem, & maxnpe, nelem, node, np, npe, xc, yc ) end if end do end do end if end do end if end do return end subroutine cross ( px, py, qx, qy, sl, sm, sh, scon, xl, xm, xh, yl, ym, yh ) !*****************************************************************************80 ! !! CROSS finds the two places where the value S = SCON occurs on a triangle. ! ! Discussion: ! ! The corners of the triangle are (XL,YL), (XM,YM) and ! (XH,YH), and the associated S values are SL, SM and SH. It ! must be the case that SL < = SM <= SH. ! ! CROSS returns two points: ! ! (PX,PY), which occurs on one of the two sides that include (XM,YM), and ! (QX,QY), which occurs on the side between (XL,YL) and (XH,YH). ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real PX, PY, the X and Y coordinates of a point ! at which S = SCON, lying on a side of the triangle which ! ends at (XM,YM). ! ! Output, real QX, QY, the X and Y coordinates of a point ! at which S = SCON, lying on the side of the triangle which ! lies between (XL,YL) and (XH,YH). ! ! Input, real SL, SM, SH, the low, medium, and high values ! of S, associated with the three corners. ! ! Input, real SCON, the value of S for which a contour line is sought. ! ! Input, real XL, XM, XH, the X coordinates of the nodes ! at which the low, medium and high values of S occur. ! ! Input, real YL, YM, YH, the Y coordinates of the nodes ! at which the low, medium and high values of S occur. ! implicit none real px real py real qx real qy real scon real sl real sm real sh real xl real xm real xh real yl real ym real yh if ( scon < sl ) then px = xl py = yl qx = xl qy = yl else if ( scon >= sh ) then px = xh py = yh qx = xh qy = yh else if ( scon < sm ) then px = xl + ( scon - sl ) * ( xm - xl ) / ( sm - sl ) py = yl + ( scon - sl ) * ( ym - yl ) / ( sm - sl ) else px = xm + ( scon - sm ) * ( xh - xm ) / ( sh - sm ) py = ym + ( scon - sm ) * ( yh - ym ) / ( sh - sm ) end if qx = xl + ( scon - sl ) * ( xh - xl ) / ( sh - sl ) qy = yl + ( scon - sl ) * ( yh - yl ) / ( sh - sl ) end if return end subroutine cross_l2 ( icross, s1, s2, scon, x1, x2, xc, y1, y2, yc ) !*****************************************************************************80 ! !! CROSS_L2 seeks contour crossings along a line. ! ! Discussion: ! ! CROSS_L2 is given the coordinates of two points in the plane, and ! the values of a quantity S there. It is assumed that S varies ! linearly along the line through the two points. ! ! CROSS_L2 is also given a desired value SC of S, and is supposed ! to seek a point, lying between the two given points, at which ! this value of S is achieved. ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ICROSS. ! ! -3, X1 = X2 and Y1=Y2. ! -2, S1 and S2 and SC are equal. ! -1, S1 and S2 are equal, and no crossing was found. ! ! 0, S1 and S2 are distinct, and no crossing was found. ! ! 1, S1 and S2 are distinct, and a crossing was found. ! 2, S1 and S2 are distinct, and S1 = SC. ! 3, S2 and S2 are distinct, and S2 = SC. ! ! Input, real S1, S2, the values of S at (X1,Y1) and (X2,Y2). ! ! Input, real SCON, the S value at which a crossing is sought. ! ! Input, real X1, X2, the X coordinates of the two endpoints. ! ! Output, real XC. ! If ICROSS = 1, 2 or 3, then XC is the X coordinate of the crossing. ! Otherwise XC is 0. ! ! Input, real Y1, Y2, the Y coordinates of the two endpoints. ! ! Output, real YC. ! If ICROSS = 1, 2, or 3, then YC is the Y coordinate of the crossing. ! Otherwise YC is 0. ! implicit none integer icross logical inside real s1 real s2 real scon real x1 real x2 real xc real y1 real y2 real yc xc = 0.0E+00 yc = 0.0E+00 ! ! Are the two points distinct? ! if ( x1 == x2 .and. y1 == y2 ) then icross = -3 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CROSS_L2 - Fatal error!' write ( *, * ) ' X1 = X2 = ', x1 write ( *, * ) ' Y1 = Y2 = ', y1 stop end if ! ! Are the two S values distinct? ! if ( s1 == s2 ) then if ( scon == s1 ) then icross = -2 else icross = -1 end if return end if ! ! Where does SCON fit in? ! call line_seg_contains_point_1d ( s1, scon, s2, inside ) if ( s2 == scon ) then icross = 3 xc = x2 yc = y2 else if ( s1 == scon ) then icross = 2 xc = x1 yc = y1 else if ( inside ) then icross = 1 xc = x1 + ( scon - s1 ) * ( x2 - x1 ) / ( s2 - s1 ) yc = y1 + ( scon - s1 ) * ( y2 - y1 ) / ( s2 - s1 ) else icross = 0 end if return end subroutine cross_l3 ( s1, s2, s3, scon, x1, x2, x3, y1, y2, y3 ) !*****************************************************************************80 ! !! CROSS_L3 seeks contour crossings in a 3 node linear element. ! ! Discussion: ! ! The coordinates of the corners of a triangle and the values of a ! scalar quantity S at those corners are given, as well as SCON, a ! special value of S. ! ! Linear interpolation is used to extend the definition of S ! to the entire triangle, and then a contour line is drawn ! to represent those locations within the triangle which ! would have the value SCON. ! ! Diagram: ! ! The triangular element has the logical form: ! ! 2 ! |\ ! | \ ! | \ ! 3---1 ! ! Modified: ! ! 07 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real S1, S2, S3, the values of S at the corners. ! ! Input, real SCON, the contour value of S. ! ! Input, real X1, X2, X3, the values of X at the corners. ! ! Input, real Y1, Y2, Y3, the values of Y at the corners. ! implicit none integer icross integer npts real s1 real s2 real s3 real scon real x1 real x2 real x3 real xc real xpts(3) real y1 real y2 real y3 real yc real