PROGRAM NSGFILE
C=======================================================================
C FINITE ELEMENT GRAPHICS FOR TWO-DIMENSIONAL CFD SOLUTIONS
C ELEMENT TYPE: FOUR-NODED ISOPARAMETRIC ELEMENT
C CREATE FILES CONTAINING XY-COORDINTE VALUES
C ORIGINAL CODE: NSEQ8GRA.FOR, EIJI FUKUMORI, JULY 1985
C=======================================================================
PARAMETER ( NF=9, MXE=90000,MXN=99000,MXB=30000,MXW=110,
* ND=4,ND2=ND*ND,INTEPT=2, INPT0=1 )
IMPLICIT REAL*8 ( A-H , O-Z )
C ARRAYS
CHARACTER*11 FNAME(NF), NAME
DIMENSION SK(ND,ND),X(ND),Y(ND),BX(2,ND), S(ND)
DIMENSION SAI(INTEPT), W(INTEPT), BPP(2,ND,INTEPT,INTEPT),
* SF(ND,INTEPT,INTEPT)
DIMENSION SAI0(INPT0), W0(INPT0), BP0(2,ND,INPT0,INPT0)
DIMENSION AM(MXN,MXW)
DIMENSION NODEX(MXE,ND),ISEG(MXE,ND),IB(MXE,ND)
DIMENSION IBNDS(MXB),BVS(MXB),IBNDFX(MXB),IBNDFY(MXB),IBNDT(MXB),
* BVX(MXB),BVY(MXB),BVT(MXB)
DIMENSION XCOORD(MXN),YCOORD(MXN),U(MXN),V(MXN),T(MXN),STM(MXN),
* P(MXN)
COMMON / DOMAIN / XMIN, XMAX, YMIN , YMAX
C=======================================================================
CALL GRULE ( INTEPT, SAI, W )
CALL DERIV ( ND, INTEPT, X, Y, SAI, BPP )
CALL SHAPEF( ND, INTEPT, X, SAI, SF )
CALL GRULE ( INPT0, SAI0, W0 )
CALL DERIV ( ND, INPT0, X, Y, SAI0, BP0 )
WRITE(*,*)' CFD GRAPHICS PROGRAM'
WRITE (*,*)' READING IN DATA FILES'
CALL INPUT ( MXE,MXN,MXB,ND,NE,NNODE,NBFX,NBFY,NBT,VISCO,
* FLMDA,NODEX,XCOORD,YCOORD,U,V,T,IBNDFX,IBNDFY,IBNDT,BVX,BVY,BVT,
* NBS,IBNDS,BVS,NF,FNAME,NAME,TIME,TREF )
WRITE (*,*)' PROJECT NAME =======>', NAME
CALL MINMAX ( MXN, NNODE, XCOORD, XMIN , XMAX )
CALL MINMAX ( MXN, NNODE, YCOORD, YMIN , YMAX )
WRITE (*,*)' MAX & MIN IN X-COORDINATE: ', XMAX, XMIN
WRITE (*,*)' MAX & MIN IN Y-COORDINATE: ', YMAX, YMIN
IF ( XMAX-XMIN .LE. 0. ) STOP' OBJECT HAS NO LENGTH IN XCOORD.'
IF ( YMAX-YMIN .LE. 0. ) STOP' OBJECT HAS NO LENGTH IN YCOORD.'
C=======================================================================
CALL BANDWID ( MXE, ND, NE, NODEX, NBW )
IF ( NBW .GT. MXW ) STOP' ERROR #5'
CALL COMPP ( MXE,MXN,INPT0,ND,BP0,NE, NNODE,FLMDA,BX,XCOORD,
* YCOORD,U,V,NODEX,IB,P)
CALL COMPSTM ( MXE,MXN,MXW,INTEPT,ND,BPP,W,SF,NE,NNODE,NBW,XCOORD,
* YCOORD,NODEX,SK,BX,U,V,AM,STM,MXB,NBS,BVS,IBNDS)
CALL CHKMAX ( MXN, NNODE, U, V, P )
NMENU = 7
ITYPE = 2
10 CALL MENU ( NBT, ITYPE )
IF ( ITYPE .EQ. 1 ) CALL PLTEL4 ( ND,MXE,MXN,NODEX,NE,
* XCOORD,YCOORD, TIME, NNODE,IB )
IF ( ITYPE .EQ. 2 ) CALL PLTUV ( MXN,NNODE,XCOORD,YCOORD,U,V,NE,
* MXE,ND,NODEX,IB,TIME,VISCO,FLMDA )
IF ( ITYPE .EQ. 3 ) CALL STREAM ( MXE,MXN,ND,NE,NNODE,NODEX,
* XCOORD,YCOORD,S,BX,IB,STM,TIME )
IF ( ITYPE .EQ. 4 ) CALL PRESS (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,
* YCOORD, P,S,BX, IB, TIME )
IF ( ITYPE .EQ. 5 ) CALL PRESSX (MXE,MXN,MXB,ND,NE,NNODE,NODEX,
* XCOORD,YCOORD, P,S,BX, IB, TIME,NBFX,IBNDFX )
IF ( ITYPE .EQ. 6 ) CALL PLTDISP ( MXN,NNODE,XCOORD,YCOORD,U,V,NE,
* MXE,ND,NODEX,IB,TIME,VISCO,FLMDA )
IF ( ITYPE .EQ. 7 ) CALL TEMP ( MXE,MXN,ND,NE,NNODE,NODEX,
* XCOORD,YCOORD,T,S,BX,IB,TREF,TIME )
IF (ITYPE.GT.NMENU ) STOP 'TERMINATION'
IF (ITYPE.LE.0) STOP 'TERMINATION'
GO TO 10
END
C
C
SUBROUTINE MENU ( NBT, ITYPE )
WRITE (*,*)'----------------------------------------------------'
WRITE (*,*)' ID # MENU'
WRITE (*,*)'===================================================='
WRITE (*,*)' 0 OR LESS TERMINATION'
WRITE (*,*)' 1 FINITE ELEMENT DISCRETIZATION'
WRITE (*,*)' 2 VELOCITY VECTOR'
WRITE (*,*)' 3 STREAM FUNCTION'
WRITE (*,*)' 4 PRESSURE CONTOUR'
WRITE (*,*)' 5 PRESSURE IN X-COORDINATE'
WRITE (*,*)' 6 DISPLACEMENT FIELD'
IF ( NBT .NE. 0 ) THEN
WRITE (*,*)' 7 TEMPERATURE CONTOUR'
ENDIF
WRITE (*,*)' 8 OR MORE TERMINATION'
WRITE (*,*)'----------------------------------------------------'
WRITE (*,*)' TYPE IN ID # = '
READ (*,*) ITYPE
RETURN
END
C
C
SUBROUTINE GRULE ( N, SAI, W )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(N), W(N)
IF ( N .LT. 1 ) STOP'N<1'
IF ( N .GT. 3 ) STOP'N>3'
IF ( N .EQ. 1 ) THEN
SAI(1) = 0.
W(1) = 2.
ENDIF
IF ( N .EQ. 2 ) THEN
SAI(1) = 1./ DSQRT(3.D0)
W(1) = 1.
SAI(2) = -SAI(1)
W(2) = W(1)
ENDIF
IF ( N .EQ. 3 ) THEN
SAI(1) = DSQRT(3.D0/5.D0)
W(1) = 5./ 9.
SAI(2) = 0.
W(2) = 8./ 9.
SAI(3) = -SAI(1)
W(3) = W(1)
ENDIF
RETURN
END
C
C
SUBROUTINE DERIV ( ND, INTEPT, BP, F, SAI, BPP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION SAI(INTEPT),BPP(2,ND,INTEPT,INTEPT), BP(ND),F(ND)
C------- SAI AND W ARE COORDINATES OF INTEGRATION POINTS AND
C WEIGHTING FACTORS.
C
DO 40 K = 1 , INTEPT
E1 = SAI(K)
DO 30 L = 1 , INTEPT
E2 = SAI(L)
C------- COMPUTATION OF BP(J) = D N(J) / D ETA1
X = E1 + 0.5
CALL ISOPARA ( ND , X , E2 , F )
X = E1 - 0.5
CALL ISOPARA ( ND , X , E2 , BP )
DO 10 I = 1 , ND
BPP(1,I,K,L) = F(I) - BP(I)
10 CONTINUE
C------- COMPUTATION OF BP(J) = D N(J) / D ETA2
Y = E2 + 0.5
CALL ISOPARA ( ND , E1 , Y , F )
Y = E2 - 0.5
CALL ISOPARA ( ND , E1 , Y , BP )
DO 20 I = 1 , ND
BPP(2,I,K,L) = F(I) - BP(I)
20 CONTINUE
C
30 CONTINUE
40 CONTINUE
RETURN
END
C
C
SUBROUTINE SHAPEF ( ND , INTEPT , F , SAI , SF )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND) , SAI(INTEPT) , SF(ND,INTEPT,INTEPT)
DO 40 K = 1 , INTEPT
E1 = SAI(K)
DO 30 L = 1 , INTEPT
E2 = SAI(L)
CALL ISOPARA ( ND , E1 , E2 , F )
DO 20 I = 1 , ND
SF(I,K,L) = F(I)
20 CONTINUE
30 CONTINUE
40 CONTINUE
RETURN
END
C
C
SUBROUTINE ISOPARA ( ND , E1 , E2 , F )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION F(ND)
F(1) = 0.25 * (1.-E1)*(1.-E2)
F(2) = 0.25 * (1.+E1)*(1.-E2)
F(3) = 0.25 * (1.+E1)*(1.+E2)
F(4) = 0.25 * (1.-E1)*(1.+E2)
RETURN
END
C
C
SUBROUTINE FORM ( MXN,MXB,MXW,NNODE,NB,NBW,A,RHS,BV,IBND )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION RHS(MXN),A(MXN,MXW),BV(MXB),IBND(MXB)
DO 50 K = 1 , NB
I = IBND(K)
DO 20 J = 2 , NBW
I = I - 1
IF ( I .LE. 0 ) GO TO 30
RHS(I) = RHS(I) - BV(K) * A(I,J)
20 CONTINUE
30 I = IBND(K)
DO 40 J = 2 , NBW
I = I + 1
IF ( I .GT. NNODE ) GO TO 50
RHS(I) = RHS(I) - BV(K) * A(IBND(K),J)
40 CONTINUE
50 CONTINUE
DO 70 K = 1 , NB
I = IBND(K)
RHS(I) = BV(K)
A(I,1) = 1.
DO 60 J = 2 , NBW
L = I - J + 1
A(I,J) = 0.
IF ( L .LE. 0 ) GO TO 60
A(L,J) = 0.
60 CONTINUE
70 CONTINUE
RETURN
END
C
C
SUBROUTINE SYSTEMQ ( MXN , MXW , NNODE , NBW , A , RHS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,MXW), RHS(MXN)
DO 30 N = 1 , NNODE
DO 20 L = 2 , NBW
C = A(N,L) / A(N,1)
I = N + L - 1
IF ( I .GT. NNODE ) GO TO 20
J = 0
DO 10 K = L , NBW
J = J + 1
A(I,J) = A(I,J) - C*A(N,K)
10 CONTINUE
A(N,L) = C
RHS(I) = RHS(I) - A(N,L) * RHS(N)
20 CONTINUE
RHS(N) = RHS(N) / A(N,1)
30 CONTINUE
C-------- BACK SUBSTITUTION
40 DO 50 K = 2 , NBW
L = N + K - 1
IF ( L .GT. NNODE ) GO TO 60
RHS(N) = RHS(N) - A(N,K)*RHS(L)
50 CONTINUE
60 N = N - 1
IF ( N .GT. 0 ) GO TO 40
RETURN
END
C
C
SUBROUTINE INPUT ( MXE,MXN,MXB,ND,NE,NNODE,NBFX,NBFY,NBT,VISCO,
* FLMDA,NODEX,XCOORD,YCOORD,U,V,T,IBNDFX,IBNDFY,IBNDT,BVX,BVY,BVT,
* NBS,IBNDS,BVS,NF,FNAME,NAME,TIME,TREF )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBNDS(MXB),U(MXN),
* V(MXN),T(MXN),IBNDFX(MXB),IBNDFY(MXB),BVX(MXB),BVY(MXB),BVS(MXB),
* IBNDT(MXB),BVT(MXB)
CHARACTER FNAME(NF)*11, NAME*11, FNAME8*11
COMMON / XVELO / UMIN, UMAX /YVELO/ VMIN, VMAX /FILE8/ FNAME8
COMMON / PRESSU/ PMIN, PMAX /STRM/ SMIN,SMAX
LOGICAL YES
C========> FILENAME NSDATA.FLN
INQUIRE ( FILE='NSDATA.FLN', EXIST=YES )
IF ( .NOT.YES ) STOP 'FILE DOES NOT EXIST'
OPEN ( 1, FILE='NSDATA.FLN', STATUS= 'OLD' )
READ(1,*)
READ(1,'(A11)') NAME
DO I = 1 , NF
READ(1,'(A11)') FNAME(I)
END DO
FNAME8 = FNAME(8)
CLOSE (1)
C========> FILENAME XXXXXXX.JNK
INQUIRE ( FILE=FNAME(1), EXIST=YES )
IF ( .NOT.YES ) STOP 'FILE DOES NOT EXIST'
OPEN ( 1, FILE=FNAME(1), STATUS='OLD' )
READ (1,'(4G15.7)') VISCO, FLMDA,TREF
READ (1,'(2I10)' ) NE, NNODE
READ (1,'(3I10)' ) NBFX, NBFY, NBT
READ (1,'(3G15.7)') DT, DTMAX, TMAX
CLOSE (1)
WRITE(*,*)' INITIAL DT =', DT
IF ( NE .GT. MXE) STOP'ERROR #1'
IF ( NNODE .GT. MXN) STOP'ERROR #2'
IF ( NBFX .GT. MXB) STOP'ERROR #3'
IF ( NBFY .GT. MXB) STOP'ERROR #4'
IF ( NBT .GT. MXB) STOP'ERROR #6'
IF ( NE .LE. 0 ) STOP'NE=0'
IF ( NNODE .LE. 0 ) STOP'NNODE=0'
IF ( NBFX .LE. 0 ) STOP'NBFX=0'
IF ( NBFY .LE. 0 ) STOP'NBFY=0'
IF ( NBT .LE. 0 ) WRITE(*,*)' SOL OF NSEQ'
C========> FILENAME XXXXXXX.ELE
INQUIRE ( FILE=FNAME(2), EXIST=YES )
IF ( .NOT.YES ) STOP 'FILE DOES NOT EXIST'
OPEN ( 1, FILE=FNAME(2), STATUS='OLD' )
DO I = 1 , NE
READ (1,*) IEL, ( NODEX(IEL,J), J = 1 , ND )
END DO
CLOSE (1)
C========> FILENAME XXXXXXX.NOD
INQUIRE ( FILE=FNAME(3), EXIST=YES )
IF ( .NOT.YES ) STOP 'FILE DOES NOT EXIST'
OPEN ( 1, FILE=FNAME(3), STATUS='OLD' )
DO I = 1 , NNODE
READ (1,*) NODE, XCOORD(NODE) , YCOORD(NODE)
END DO
CLOSE (1)
C========> FILENAME XXXXXXX.BND
INQUIRE ( FILE=FNAME(4), EXIST=YES )
IF ( .NOT.YES ) STOP 'FILE DOES NOT EXIST'
OPEN ( 1, FILE=FNAME(4), STATUS='OLD' )
DO I = 1 , NBFX
READ (1,*) IBNDFX(I) , BVX(I)
END DO
DO I = 1 , NBFY
READ (1,*) IBNDFY(I) , BVY(I)
END DO
IF ( NBT .NE. 0 ) THEN
DO I = 1 , NBT
READ (1,*) IBNDT(I) , BVT(I)
END DO
ENDIF
CLOSE (1)
C========> FILENAME XXXXXXX.BIN
OPEN ( 1, FILE=FNAME(6), STATUS='OLD',FORM='UNFORMATTED' )
READ (1) TIME, DT
READ (1) ( U(I) , I = 1 , NNODE )
READ (1) ( V(I) , I = 1 , NNODE )
IF ( NBT .GT. 0 ) READ (1) ( T(I) , I = 1 , NNODE )
CLOSE (1)
WRITE (*,*)' DT, DTMAX =', DT, DTMAX
WRITE (*,*)' TIME, TMAX =', TIME, TMAX
C========> FILENAME XXXXXXX.MAX
INQUIRE ( FILE=FNAME(8), EXIST=YES )
IF ( .NOT.YES ) STOP 'FILE DOES NOT EXIST'
OPEN ( 1, FILE=FNAME(8), STATUS='OLD' )
READ (1,*) UMIN , UMAX
READ (1,*) VMIN , VMAX
READ (1,*) PMIN , PMAX
READ (1,*) SMIN , SMAX
CLOSE (1)
C========> FILENAME XXXXXXX.STM
INQUIRE ( FILE=FNAME(9), EXIST=YES )
IF ( .NOT.YES ) STOP 'FILE DOES NOT EXIST'
OPEN ( 1, FILE=FNAME(9), STATUS='OLD' )
READ (1,*) NBS
READ (1,*) ( IBNDS(I), BVS(I) , I = 1 , NBS )
CLOSE (1)
RETURN
END
C
C
SUBROUTINE CHKMAX ( MXN , NNODE , U , V , P )
IMPLICIT REAL*8 ( A-H , O-Z )
CHARACTER*11 FNAME8
DIMENSION U(MXN) , V(MXN) , P(MXN)
COMMON / XVELO / UMIN, UMAX /YVELO/ VMIN, VMAX /FILE8/ FNAME8
COMMON / PRESSU/ PMIN, PMAX /STRM/ SMIN,SMAX
CALL MINMAX ( MXN, NNODE, P , PPMIN , PPMAX )
IF ( PPMIN .LT. PMIN ) PMIN = PPMIN
IF ( PPMAX .GT. PMAX ) PMAX = PPMAX
CALL MINMAX ( MXN, NNODE, U , UUMIN , UUMAX )
IF ( UUMIN .LT. UMIN ) UMIN = UUMIN
IF ( UUMAX .GT. UMAX ) UMAX = UUMAX
CALL MINMAX ( MXN, NNODE, V , VVMIN , VVMAX )
IF ( VVMIN .LT. VMIN ) VMIN = VVMIN
IF ( VVMAX .GT. VMAX ) VMAX = VVMAX
OPEN ( 1 , FILE=FNAME8, STATUS = 'OLD' )
WRITE (1,*) UMIN , UMAX
WRITE (1,*) VMIN , VMAX
WRITE (1,*) PMIN , PMAX
WRITE (1,*) SMIN , SMAX
CLOSE (1)
RETURN
END
C
C
SUBROUTINE BANDWID ( MXE , ND , NE , NODEX , NBW )
DIMENSION NODEX(MXE,ND)
NBW = 0
DO I = 1 , NE
DO J = 1 , ND-1
DO K = J+1 , ND
NBW = MAX0 ( NBW , IABS(NODEX(I,J)-NODEX(I,K)) )
END DO
END DO
END DO
NBW = NBW + 1
WRITE (*,*) ' HALH BANDWIDTH =', NBW
RETURN
END
C
C
SUBROUTINE COMPSTM ( MXE,MXN,MXW,INTEPT,ND,BPP,W,SF,NE,NNODE,NBW,
* XCOORD,YCOORD,NODEX,STIFF,B,U,V,STRM,RHS,MXB,NBS,BVS,IBNDS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* SF(ND,INTEPT,INTEPT),STRM(MXN,MXW),STIFF(ND,ND),B(2,ND),U(MXN),
* V(MXN),BVS(MXB),IBNDS(MXB),RHS(MXN),BPP(2,ND,INTEPT,INTEPT),
* W(INTEPT)
C-------- RESET
VISCO = 1.
DO 10 I = 1 , NNODE
RHS(I) = 0.
DO 10 J = 1 , NBW
STRM(I,J) = 0.
10 CONTINUE
DO 500 IEL = 1 , NE
DO 33 I = 1 , ND
DO 33 J = 1 , ND
STIFF(I,J) = 0.
33 CONTINUE
C------- GAUSS INTEGRATION
DO 400 K = 1 , INTEPT
DO 300 L = 1 , INTEPT
WEIGHT = W(K) * W(L)
YAC11 = 0.
YAC12 = 0.
YAC21 = 0.
YAC22 = 0.
DO I = 1 , ND
NODE = NODEX(IEL,I)
YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODE)
YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODE)
YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODE)
YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODE)
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
Y11 = YAC22 / DETJ
Y12 = -YAC12 / DETJ
Y21 = -YAC21 / DETJ
Y22 = YAC11 / DETJ
DO J = 1 , ND
B(1,J) = Y11 * BPP(1,J,K,L) + Y12 * BPP(2,J,K,L)
B(2,J) = Y21 * BPP(1,J,K,L) + Y22 * BPP(2,J,K,L)
END DO
C-------- COMPUTATION OF ROTATION
DUDY = 0
DVDX = 0.
DO 13 J = 1 , ND
NODE = NODEX(IEL,J)
DUDY = DUDY + B(2,J) * U(NODE)
DVDX = DVDX + B(1,J) * V(NODE)
13 CONTINUE
ROTA = ( DUDY - DVDX ) * WEIGHT * DETJ
DO 12 J = 1 , ND
NODE = NODEX(IEL,J)
RHS(NODE) = RHS(NODE) - SF(J,K,L) * ROTA
12 CONTINUE
C------ LOCAL STIFFNESS MATRIX ASEMBLY
BETA = VISCO * WEIGHT * DETJ
DO 11 I = 1 , ND-1
DO 11 J = I+1 , ND
ENTRY = ( B(1,I)*B(1,J) + B(2,I)*B(2,J) ) * BETA
STIFF(I,J) = STIFF(I,J) + ENTRY
STIFF(I,I) = STIFF(I,I) - ENTRY
STIFF(J,J) = STIFF(J,J) - ENTRY
STIFF(J,I) = STIFF(I,J)
11 CONTINUE
300 CONTINUE
400 CONTINUE
C--------- GLOBAL STIFFNESS MATRIX ASSEMBLY
DO 30 K = 1 , ND
I = NODEX(IEL,K)
DO 23 L = 1 , ND
J = NODEX(IEL,L) - I + 1
IF ( J .LE. 0 ) GO TO 23
STRM(I,J) = STRM(I,J) + STIFF(K,L)
23 CONTINUE
30 CONTINUE
500 CONTINUE
C--------- STREAM FUNCTION VALUE EVALUATION
CALL FORM ( MXN, MXB,MXW,NNODE,NBS,NBW,STRM,RHS,BVS,IBNDS )
CALL SYSTEMQ ( MXN, MXW, NNODE, NBW, STRM, RHS )
RETURN
END
C
C
SUBROUTINE STREAM ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* S,B,IB,RHS,TIME )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND),
* B(2,ND),RHS(MXN)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / PORT / STMMIN, STMMAX / INCREM / DSTM
FILENAME = "STREAM"
CALL MINMAX ( MXN, NNODE, RHS, STMMIN, STMMAX )
CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,S,B,IB,RHS)
WRITE(*,200) TIME, STMMIN, STMMAX, DSTM
200 FORMAT (1X,'STREAM FUNCTION:'/1X,'BLUE =ZERO'/1X,'RED =CLKWISE'/
* 1X, 'GREEN=A-CLKWISE' ,20(/) ,1X, 'CURRENT VALUES:'/
* 1X, 'TIME=',G10.3 /
* 1X, 'MIN =',G10.3 /1X, 'MAX =',G10.3 /1X,'DS =',G10.3 )
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE TEMP ( MXE,MXN,ND,NE,NNODE,NODEX, XCOORD,YCOORD,
* T,S,B,IB,TREF,TIME )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND),
* B(2,ND),T(MXN)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / PORT / TMIN, TMAX / INCREM / DT
FILENAME = "TEMP"
DO I = 1 , NNODE
T(I) = T(I) - TREF
END DO
CALL MINMAX ( MXN, NNODE, T, TMIN, TMAX )
CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,S,B,IB,T)
DO I = 1 , NNODE
T(I) = T(I) + TREF
END DO
CALL MINMAX ( MXN, NNODE, T, TMIN, TMAX )
WRITE(*,200) TREF,TIME,TMIN, TMAX, DT
CALL PLTEXT
200 FORMAT ( 1X,'TEMPERATURE FIELD:'/ 1X,'BLUE=TREF'/
* 1X, 'RED=BELOW TREF'/ 1X, 'GREEN=ABOVE TREF'/1X,'TREF=',G10.3,
* 19(/),1X, 'CURRENT VALUES:'/1X,'TIME=',G10.3 / 1X, 'MIN=',G10.3/
* 1X, 'MAX=',G10.3 / 1X,'DS=',G10.3 )
RETURN
END
C
C
SUBROUTINE PLTSAI ( DS, NSTEP, START, CRD, SS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION CRD(2,4), SS(4), X(4), Y(4), S(4)
IF ( NSTEP .EQ. 0 ) RETURN
X(3) = ( CRD(1,1) + CRD(1,2) + CRD(1,3) + CRD(1,4) ) / 4.
Y(3) = ( CRD(2,1) + CRD(2,2) + CRD(2,3) + CRD(2,4) ) / 4.
S(3) = ( SS(1) + SS(2) + SS(3) + SS(4) ) / 4.
SMAX = DMAX1 ( SS(1), SS(2), SS(3), SS(4) )
SMIN = DMIN1 ( SS(1), SS(2), SS(3), SS(4) )
DO 52 LEVEL = 1 , NSTEP
SXY = START + (LEVEL-1) * DS
IF ( (SMAX-SXY)*(SMIN-SXY) .LT. 0 ) THEN
IF ( SXY .LT. 0. ) CALL JCOLOR ( 9 )
IF ( SXY .GT. 0. ) CALL JCOLOR ( 10 )
IF ( SXY .EQ. 0. ) CALL JCOLOR ( 12 )
DO 60 IEL = 1 , 4
X(1) = CRD(1,IEL)
Y(1) = CRD(2,IEL)
S(1) = SS(IEL)
IF ( IEL .EQ. 4 ) THEN
X(2) = CRD(1,1)
Y(2) = CRD(2,1)
S(2) = SS(1)
ELSE
X(2) = CRD(1,IEL+1)
Y(2) = CRD(2,IEL+1)
S(2) = SS (IEL+1)
ENDIF
X(4) = X(1)
Y(4) = Y(1)
S(4) = S(1)
K = 0
DO 70 ISG = 1 , 3
IF ( S(ISG ) .LT. SXY ) GO TO 30
IF ( S(ISG+1) .LT. SXY ) GO TO 40
GO TO 70
30 IF ( S(ISG+1) .LT. SXY ) GO TO 70
40 T = ( SXY - S(ISG) ) / ( S(ISG+1) - S(ISG) )
X0 = X(ISG+1)*T + (1.- T)*X(ISG)
Y0 = Y(ISG+1)*T + (1.- T)*Y(ISG)
IF ( K .EQ. 0 ) GO TO 71
CALL XDRAW ( X0, Y0 )
GO TO 60
71 CALL XMOVE ( X0, Y0 )
K = 1
70 CONTINUE
60 CONTINUE
ENDIF
52 CONTINUE
RETURN
END
C
C
SUBROUTINE PLTUV (MXN,NNODE,XCOORD,YCOORD,U,V,NE,MXE,ND,NODEX,
* IB, TIME, VISCO, FLMDA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(MXN),YCOORD(MXN),U(MXN),V(MXN),NODEX(MXE,ND),
* IB(MXE,ND)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / DOMAIN / XMIN, XMAX, YMIN, YMAX
FILENAME = "VECTOR"
PLTLMT = 0.001
SQLMT = PLTLMT **2
RATIO = 0.05
RMAX = 0.
DO I = 1 , NNODE
RMAX = DMAX1 ( RMAX, U(I)*U(I)+V(I)*V(I) )
END DO
IF ( RMAX .EQ. 0. ) RETURN
RMAX = DSQRT ( RMAX )
CALL PLTLGO
CALL MINMAX ( MXN, NNODE, U, UMIN, UMAX )
CALL MINMAX ( MXN, NNODE, V, VMIN, VMAX )
WRITE(*,200) TIME, VISCO, FLMDA, UMIN, UMAX, VMIN, VMAX
CALL JCOLOR ( 11 )
IC = 0
CALL BOUND (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,IB,IC )
CALL JCOLOR ( 11 )
FACT = RATIO * DMAX1 ( (YMAX-YMIN) , (XMAX-XMIN) ) / RMAX
DO I = NNODE , 1 , -1
IF ( (U(I)*U(I)+V(I)*V(I))/RMAX .GT. SQLMT )
* CALL VECPLT ( XCOORD(I),YCOORD(I),U(I),V(I),FACT )
END DO
CALL PLTEXT
200 FORMAT ( 1X, 'VELOCITY'/1X,'VECTOR PLOT'/
* 1X,'CURRENT VALUES:'/1X,'TIME=',G10.3/1X,'VISCO=',G10.3/
* 1X,'FLMDA=',G10.3,19(/),1X,'UMIN=',G10.3/
* 1X,'UMAX=',G10.3/ 1X,'VMIN=',G10.3/1X, 'VMAX=',G10.3 )
RETURN
END
C
C
SUBROUTINE VECPLT ( X0, Y0, U, V, FACT )
IMPLICIT REAL*8 ( A-H , O-Z )
DATA AL, BETA / 0.9D0, 0.08D0 /
CALL XMOVE ( X0, Y0 )
DX = U * FACT
DY = V * FACT
CALL XDRAW ( X0+AL*DX , Y0+AL*DY )
CALL XMOVE ( X0+AL*DX , Y0+AL*DY )
RNX = BETA * FACT * V
RNY = - BETA * FACT * U
CALL XDRAW ( X0+AL*DX+RNX , Y0+AL*DY+RNY )
CALL XMOVE ( X0+AL*DX+RNX , Y0+AL*DY+RNY )
CALL XDRAW ( X0+ DX , Y0+ DY )
CALL XMOVE ( X0+ DX , Y0+ DY )
CALL XDRAW ( X0+AL*DX-RNX , Y0+AL*DY-RNY )
CALL XMOVE ( X0+AL*DX-RNX , Y0+AL*DY-RNY )
CALL XDRAW ( X0+AL*DX , Y0+AL*DY )
RETURN
END
C
C
SUBROUTINE PLTDISP (MXN,NNODE,XCOORD,YCOORD,U,V,NE,MXE,ND,NODEX,
* IB, TIME, VISCO, FLMDA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(MXN),YCOORD(MXN),U(MXN),V(MXN),NODEX(MXE,ND),
* IB(MXE,ND)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / DOMAIN / XMIN, XMAX, YMIN, YMAX
DATA AL, BETA / 0.9D0, 0.08D0 /
FILENAME = "DISPLACE"
PLTLMT = 0.001
SQLMT = PLTLMT **2
RATIO = 0.05
RMAX = 0.
DO I = 1 , NNODE
RMAX = DMAX1 ( RMAX, U(I)*U(I)+V(I)*V(I) )
END DO
IF ( RMAX .EQ. 0. ) RETURN
RMAX = DSQRT ( RMAX )
CALL PLTLGO
CALL MINMAX ( MXN, NNODE, U, UMIN, UMAX )
CALL MINMAX ( MXN, NNODE, V, VMIN, VMAX )
FACT = RATIO * DMAX1 ( (YMAX-YMIN) , (XMAX-XMIN) ) / RMAX
DO IEL = 1 , NE
C------ 1ST SEGMENT
X1 = XCOORD(NODEX(IEL,1))+AL*U(NODEX(IEL,1))*FACT
Y1 = YCOORD(NODEX(IEL,1))+AL*V(NODEX(IEL,1))*FACT
X2 = XCOORD(NODEX(IEL,2))+AL*U(NODEX(IEL,2))*FACT
Y2 = YCOORD(NODEX(IEL,2))+AL*V(NODEX(IEL,2))*FACT
CALL XMOVE ( X1,Y1 )
CALL XDRAW ( X2,Y2 )
C------ 2ND SEGMENT
X1 = XCOORD(NODEX(IEL,2))+AL*U(NODEX(IEL,2))*FACT
Y1 = YCOORD(NODEX(IEL,2))+AL*V(NODEX(IEL,2))*FACT
X2 = XCOORD(NODEX(IEL,3))+AL*U(NODEX(IEL,3))*FACT
Y2 = YCOORD(NODEX(IEL,3))+AL*V(NODEX(IEL,3))*FACT
CALL XMOVE ( X1,Y1 )
CALL XDRAW ( X2,Y2 )
C------ 3RD SEGMENT
X1 = XCOORD(NODEX(IEL,3))+AL*U(NODEX(IEL,3))*FACT
Y1 = YCOORD(NODEX(IEL,3))+AL*V(NODEX(IEL,3))*FACT
X2 = XCOORD(NODEX(IEL,4))+AL*U(NODEX(IEL,4))*FACT
Y2 = YCOORD(NODEX(IEL,4))+AL*V(NODEX(IEL,4))*FACT
CALL XMOVE ( X1,Y1 )
CALL XDRAW ( X2,Y2 )
C------ 4TH SEGMENT
X1 = XCOORD(NODEX(IEL,4))+AL*U(NODEX(IEL,4))*FACT
Y1 = YCOORD(NODEX(IEL,4))+AL*V(NODEX(IEL,4))*FACT
X2 = XCOORD(NODEX(IEL,1))+AL*U(NODEX(IEL,1))*FACT
Y2 = YCOORD(NODEX(IEL,1))+AL*V(NODEX(IEL,1))*FACT
CALL XMOVE ( X1,Y1 )
CALL XDRAW ( X2,Y2 )
END DO
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE BOUND (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,IB,IC)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(MXN),YCOORD(MXN),NODEX(MXE,ND),IB(MXE,ND)
C---------- IF IC=0, BOUNDARY LINE ONLY
C---------- IF IC=1, ELEMENTS
C--------- INITIALIZATION
DO I = 1 , NE
DO J = 1 , ND
IB(I,J) = 0
END DO
END DO
C
DO IEL = 1 , NE-1
DO ISEG = 1 , ND
IF ( IB(IEL,ISEG) .EQ. 0 ) THEN
IS = ISEG
IE = ISEG+1
IF ( IS .EQ. ND ) IE = 1
IS = NODEX(IEL,IS)
IE = NODEX(IEL,IE)
DO JEL = IEL+1 , NE
DO JSEG = 1 , ND
IF ( IB(JEL,JSEG) .EQ. 0 ) THEN
JS = JSEG
JE = JSEG+1
IF ( JS .EQ. ND ) JE = 1
JS = NODEX(JEL,JS)
JE = NODEX(JEL,JE)
IF ( (IS .EQ. JE) .AND. (IE .EQ. JS) ) THEN
IF ( IC .EQ. 0 ) IB(IEL,ISEG) = 1
IB(JEL,JSEG) = 1
END IF
END IF
END DO
END DO
END IF
END DO
END DO
C
DO IEL = 1 , NE
DO ISEG = 1 , ND
IF ( IB(IEL,ISEG) .EQ. 0 ) THEN
IS = ISEG
IE = ISEG+1
IF ( IS .EQ. ND ) IE = 1
IS = NODEX(IEL,IS)
IE = NODEX(IEL,IE)
CALL XMOVE ( XCOORD(IS) , YCOORD(IS) )
CALL XDRAW ( XCOORD(IE) , YCOORD(IE) )
END IF
END DO
END DO
RETURN
END
C
C
SUBROUTINE PRESS (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* P,S,B, IB, TIME )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* P(MXN), S(ND),B(2,ND),IB(MXE,ND)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / DOMAIN / XMIN, XMAX, YMIN, YMAX / PRESSU / PMIN, PMAX
COMMON / PORT / QMIN, QMAX / INCREM / DP
FILENAME = "PRESSURE"
CALL MINMAX ( MXN, NNODE, P, QMIN, QMAX )
CALL CONTOUR ( MXE, MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,S,B,IB,P)
WRITE(*,200)
WRITE(*,210) TIME, QMIN, QMAX, DP
CALL PLTEXT
200 FORMAT ( 1X,'PRESSURE CONTOUR'/1X,'BLUE=ZERO'/1X,'RED=-'/
* 1X,'GREEN=+',19(/) )
210 FORMAT ( 1X,'CURRENT VALUES:'/1X,'TIME=',G10.3/1X,'PMIN=', G10.3/
* 1X, 'PMAX=',G10.3 / 1X,'DP=',G10.3 )
RETURN
END
C
C
SUBROUTINE PRESSX (MXE,MXN,MXB,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* P,S,B, IB, TIME,NBFX,IBNDFX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IBNDFX(MXB),
* P(MXN), S(ND),B(2,ND),IB(MXE,ND)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / DOMAIN / XMIN, XMAX, YMIN, YMAX / PRESSU / PMIN, PMAX
COMMON / PORT / QMIN, QMAX / INCREM / DP
FILENAME = "PRESSINX"
CALL MINMAX ( MXN, NNODE, P, QMIN, QMAX )
WRITE(*,200)
WRITE(*,210) TIME, QMIN, QMAX, DP
200 FORMAT ( 1X,'PRESSURE CONTOUR IN X' )
210 FORMAT ( 1X,'CURRENT VALUES:'/1X,'TIME=',G10.3/1X,'PMIN=', G10.3/
* 1X, 'PMAX=',G10.3 / 1X,'DP=',G10.3 )
C
OPEN ( 1, FILE=FILENAME, STATUS='UNKNOWN' )
DO J = 1, NBFX
I = IBNDFX(J)
WRITE(1,300) XCOORD(I), P(I)
300 FORMAT ( 2F26.14 )
END DO
CLOSE (1)
RETURN
END
C
C
SUBROUTINE CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* S, B, IB, P )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), XCOORD(MXN),YCOORD(MXN), P(MXN),
* S(ND), B(2,ND), IB(MXN)
COMMON / PORT / PPMIN, PPMAX, / INCREM / DS
IF ( PPMAX .EQ. PPMIN ) RETURN
WRITE(*,*)' TYPE IN NUMBER OF CONTOUR LINES(N=20 ABOUT RIGHT)'
WRITE(*,*)' N=0 FOR QUIT, MAX N = 99'
WRITE(*,*)' N='
READ (*,*) NSTEP
IF ( NSTEP .LE. 0 ) RETURN
IF ( NSTEP .GT. 99 ) RETURN
NSTEP = NSTEP/2*2 + 1
DS = ( PPMAX - PPMIN ) / NSTEP
C
IF ( (PPMIN .LT. 0.D0) .AND. (PPMAX.GT.0.D0) ) THEN
DS = DMAX1(PPMAX, -PPMIN)/NSTEP
END IF
C
CALL PLTLGO
CALL JCOLOR ( 9 )
IC = 0
CALL BOUND ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,IB,IC )
DO IEL = 1 , NE
DO I = 1 , ND
B(1,I) = XCOORD(NODEX(IEL,I))
B(2,I) = YCOORD(NODEX(IEL,I))
S(I) = P(NODEX(IEL,I))
END DO
IF ( (PPMIN .LT. 0.D0) .AND. (PPMAX.GT.0.D0) ) THEN
CALL PLTSAI ( DS, NSTEP, -NSTEP*DS, B, S )
CALL PLTSAI ( DS, 1, 0.D0, B, S )
CALL PLTSAI ( DS, NSTEP, DS, B, S )
ELSE
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
END IF
END DO
RETURN
END
C
C
SUBROUTINE MINMAX ( MXN, NNODE, Q, QMIN, QMAX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION Q(MXN)
QMIN = Q(1)
QMAX = Q(1)
DO I = 1 , NNODE
QMIN = DMIN1 ( QMIN , Q(I) )
QMAX = DMAX1 ( QMAX , Q(I) )
END DO
RETURN
END
C
SUBROUTINE PLTEL4 (ND,MXE,MXN,NODEX,NE,XCOORD,YCOORD,TIME,
* NNODE,IB)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),IB(MXE,ND)
CHARACTER FILENAME*8
COMMON / PLOTTING / FILENAME
COMMON / DOMAIN / XMIN, XMAX, YMIN, YMAX
FILENAME = "ELEMENT"
CALL PLTLGO
WRITE (*,200) TIME, NE, NNODE, XMIN, XMAX, YMIN, YMAX
IC = 1
CALL BOUND ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,IB,IC )
CALL PLTEXT
200 FORMAT ( 1X,'FINITE ELEMENT'/ 1X,'DISCRETIZATION'/
* 1X,'TIME=',G10.3 /1X,'NE=', I5 / 1X,'NNODE=', I5, 20(/),
* 1X,'XMIN=',G10.3 / 1X,'XMAX=',G10.3 /
* 1X,'YMIN=',G10.3 / 1X,'YMAX=',G10.3 )
RETURN
END
C
C
SUBROUTINE SEGCHK ( ND,MXE,NODEX,NE,ISEG,II,JJ )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND), ISEG(MXE,ND)
DO IEL = 1 , NE
DO I = 1 , ND
IF ( NODEX(IEL,I) .EQ.JJ ) THEN
J = I + 1
IF ( I .EQ. ND ) J = 1
IF ( NODEX(IEL,J).EQ.II ) THEN
ISEG (IEL,I) = 0
RETURN
ENDIF
ENDIF
END DO
END DO
RETURN
END
C
C
SUBROUTINE COMPP ( MXE,MXN,INTEPT,ND,BPP,NE,NNODE,FLMDA,
* BX,XCOORD,YCOORD,U,V,NODEX, IRPT, P )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN), BX(2,ND),
* BPP(2,ND,INTEPT,INTEPT),U(MXN),V(MXN), IRPT(MXN),P(MXN)
C-------- COMPUTATION OF PRESSURE
C-------- RETURN VALUES P(MXN)
C-------- RESET
DO I = 1 , NNODE
P(I) = 0.
IRPT(I) = 0
END DO
C
NUMBERP = INTEPT*INTEPT
DO 500 IEL = 1 , NE
CONTI = 0.
C------- INTEGRATION BY GAUSS
DO 400 K = 1 , INTEPT
DO 300 L = 1 , INTEPT
YAC11 = 0.
YAC12 = 0.
YAC21 = 0.
YAC22 = 0.
DO I = 1 , ND
NODE = NODEX(IEL,I)
YAC11 = YAC11 + BPP(1,I,K,L) * XCOORD(NODE)
YAC12 = YAC12 + BPP(1,I,K,L) * YCOORD(NODE)
YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(NODE)
YAC22 = YAC22 + BPP(2,I,K,L) * YCOORD(NODE)
END DO
DETJ = YAC11 * YAC22 - YAC12 * YAC21
Y11 = YAC22 / DETJ
Y12 = -YAC12 / DETJ
Y21 = -YAC21 / DETJ
Y22 = YAC11 / DETJ
DO J = 1 , ND
BX(1,J) = Y11 * BPP(1,J,K,L) + Y12 * BPP(2,J,K,L)
BX(2,J) = Y21 * BPP(1,J,K,L) + Y22 * BPP(2,J,K,L)
END DO
DO J = 1 , ND
NODE = NODEX(IEL,J)
CONTI = CONTI + (BX(1,J)*U(NODE)+BX(2,J)*V(NODE))
END DO
300 CONTINUE
400 CONTINUE
PRESS = - FLMDA*CONTI / NUMBERP
C--------- DISTRIBUTION OF PRESSURE
DO J = 1 , ND
I = NODEX(IEL,J)
P(I) = P(I) + PRESS
IRPT(I) = IRPT(I) + 1
END DO
500 CONTINUE
DO I = 1 , NNODE
P(I) = P(I) / IRPT(I)
END DO
RETURN
END
C======================== GRAPHICS RUTINES =============================
SUBROUTINE PLTLGO
IMPLICIT REAL*8 ( A-H , O-Z )
CHARACTER FILENAME*8
COMMON / DOMAIN / XMIN , XMAX , YMIN , YMAX
COMMON / PENDOMAIN / IXMIN , IXMAX , IYMIN , IYMAX
COMMON / PLOTTING / FILENAME, RMAGNI, DIGITS
DIGITS = 99999.D0
OPEN ( 1, FILE=FILENAME, STATUS='UNKNOWN' )
RMAGNI = DMAX1 ( XMAX-XMIN, YMAX-YMIN )
IXMIN = DIGITS
IYMIN = DIGITS
IXMAX = -DIGITS
IYMAX = -DIGITS
RETURN
END
C
C
SUBROUTINE PLTEXT
COMMON / PENDOMAIN / IXMIN , IXMAX , IYMIN , IYMAX
WRITE(1,'(4I7)') IXMIN , IXMAX , IYMIN , IYMAX
CLOSE (1)
RETURN
END
C
C
SUBROUTINE JCOLOR ( I )
RETURN
END
C
C
SUBROUTINE XMOVE ( X , Y )
IMPLICIT REAL*8 ( A-H , O-Z )
COMMON / DOMAIN / XMIN , XMAX , YMIN , YMAX
COMMON / PENDOMAIN / IXMIN , IXMAX , IYMIN , IYMAX
COMMON / PLOTTING / FILENAME, RMAGNI, DIGITS
COMMON / PENMOVE /IXMOVE, IYMOVE
IXMOVE = (X-XMIN)/RMAGNI*DIGITS
IYMOVE = (Y-YMIN)/RMAGNI*DIGITS
IXMIN = MIN0 ( IXMIN, IXMOVE )
IXMAX = MAX0 ( IXMAX, IXMOVE )
IYMIN = MIN0 ( IYMIN, IYMOVE )
IYMAX = MAX0 ( IYMAX, IYMOVE )
RETURN
END
C
C
SUBROUTINE XDRAW ( X , Y )
IMPLICIT REAL*8 ( A-H , O-Z )
COMMON / PENDOMAIN / IXMIN , IXMAX , IYMIN , IYMAX
COMMON / DOMAIN / XMIN , XMAX , YMIN , YMAX
COMMON / PLOTTING / FILENAME, RMAGNI, DIGITS
COMMON / PENMOVE /IXMOVE, IYMOVE
IX = (X-XMIN)/RMAGNI*DIGITS
IY = (Y-YMIN)/RMAGNI*DIGITS
WRITE(1,'(4I7)') IXMOVE, IYMOVE, IX, IY
IXMIN = MIN0 ( IXMIN, IX )
IXMAX = MAX0 ( IXMAX, IX )
IYMIN = MIN0 ( IYMIN, IY )
IYMAX = MAX0 ( IYMAX, IY )
RETURN
END