PROGRAM NSTG4FILE
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 FEB. 1, 2013
C=======================================================================
INCLUDE 'PARAM.DAT'
CCCC PARAMETER ( NF=9, MXE=1400,MXN=1500,MXB=616,MXW=38, ND=4 )
PARAMETER ( ND2=ND*ND,INPT0=INTEPT-1 )
IMPLICIT REAL*8 ( A-H , O-Z )
C ARRAYS
CHARACTER*19 FNAME(NF), NAME*15
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)
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), IRPT(MXN)
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
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,IRPT,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)
WRITE (*,*) 'AFTER CALL COMPSTM'
C=======================================================================
WRITE (*,*) 'CALL PLTEL4'
CALL PLTEL4 ( ND,MXE,MXN,NODEX,NE,XCOORD, YCOORD )
C=======================================================================
CALL PLTUV ( MXN,NNODE,XCOORD,YCOORD,U,V,NE,MXE,ND,NODEX )
CALL STREAM ( MXE,MXN,ND,NE,NNODE,NODEX,
* XCOORD,YCOORD,S,BX,STM )
CALL PRESS (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,P,S,BX )
CALL TEMP ( MXE,MXN,ND,NE,NNODE,NODEX,
* XCOORD,YCOORD,T,S,BX,TREF )
STOP 'TERMINATION'
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 K = 1 , INTEPT
E1 = SAI(K)
DO 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 I = 1 , ND
BPP(1,I,K,L) = F(I) - BP(I)
END DO
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 I = 1 , ND
BPP(2,I,K,L) = F(I) - BP(I)
END DO
C
END DO
END DO
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 K = 1 , INTEPT
E1 = SAI(K)
DO L = 1 , INTEPT
E2 = SAI(L)
CALL ISOPARA ( ND , E1 , E2 , F )
DO I = 1 , ND
SF(I,K,L) = F(I)
END DO
END DO
END DO
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 BV(MXB),IBND(MXB),RHS(MXN), A(MXN,MXW)
C---- NBWA = HALFBANDWIDTH INCLUDING DIAGONAL
C----- REFORMATION OF VECTOR {RHS} & MATRIX [A]
C------- RETURN VALUE: A, RHS
C-----REFORMATION OF VECTOR {RHS}
DO K = 1 , NB
I = IBND(K)
DO J = 2 , NBW
I = I - 1
IF ( I.GT. 0 ) RHS(I) = RHS(I) - BV(K) * A(I,J)
END DO
I = IBND(K)
DO J = 2 , NBW
I = I + 1
IF ( I .LE. NNODE ) RHS(I) = RHS(I) - BV(K) * A(IBND(K),J)
END DO
END DO
C-----REFORMATION OF MATRIX [A]
DO K = 1 , NB
I = IBND(K)
RHS(I) = BV(K)
A(I,1) = 1.
DO J = 2 , NBW
A(I,J) = 0.
IF ( I-J+1 .GT. 0 ) A(I-J+1,J) = 0.
END DO
END DO
RETURN
END
C
C
SUBROUTINE SYSTEMQ ( MXN, MXW, NUMNP, MBAND, A, B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION A(MXN,MXW) , B(MXN)
C------- RETURN VALUE: B
C---------- ELIMINATION ------------------
DO N = 1 , NUMNP
DO L = 2 , MBAND
C = A(N,L) / A(N,1)
I = N + L - 1
IF ( I .LE. NUMNP ) THEN
J = 0
DO K = L , MBAND
J = J + 1
A(I,J) = A(I,J) - C * A(N,K)
END DO
A(N,L) = C
B(I) = B(I) - A(N,L) * B(N)
ENDIF
END DO
B(N) = B(N) / A(N,1)
END DO
C---------- BACKSUBSTITUTION -------------
N = NUMNP
DO WHILE ( N .GT. 0 )
DO K = 2 , MBAND
L = N + K - 1
IF ( L .LE. NUMNP ) B(N) = B(N)-A(N,K)*B(L)
END DO
N = N - 1
END DO
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,PROJECT,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)*19, PROJECT*15
C=======================================================================
C========> FILENAME NSDATAX.FLN
IF (ND .EQ. 4) OPEN ( 1, FILE='NSDATA4.FLN', STATUS='UNKNOWN' )
IF (ND .EQ. 9) OPEN ( 1, FILE='NSDATA9.FLN', STATUS='UNKNOWN' )
READ(1,'(A15)') PROJECT
DO I = 1 , NF
READ(1,'(A19)') FNAME(I)
END DO
CLOSE (1)
WRITE(*,*)' PROJECT NAME:============>>>',PROJECT
C=======================================================================
C========> FILENAME XXXXXXX.DAT
OPEN ( 1, FILE=FNAME(1), STATUS='UNKNOWN' )
C--------> HEADER
READ (1,*) VISCO, FLMDA, TREF, CONDUCT, C1
READ (1,*) TMAX , DTMAX
READ (1,*) MAXITE, ERMAX, ITEFIX
C--------> ELEMENT CONNECTIVITY
READ (1,*) NE
IF ( NE .GT. MXE) STOP'ERROR #1'
IF ( NE .LE. 0 ) STOP'NE=0'
DO I = 1 , NE
READ (1,*) IEL, ( NODEX(IEL,J), J = 1 , ND )
END DO
C--------> NODAL INFORMATION
READ (1,*) NNODE
IF ( NNODE .GT. MXN) STOP'ERROR #2'
IF ( NNODE .LE. 0 ) STOP'NNODE=0'
DO I = 1 , NNODE
READ (1,*) NODE, XCOORD(NODE) , YCOORD(NODE)
END DO
C========> BOUNDARY CONDITIONS
READ (1,*) NBFX
DO I = 1 , NBFX
READ (1,*) IBNDFX(I) , BVX(I)
END DO
READ (1,*) NBFY
DO I = 1 , NBFY
READ (1,*) IBNDFY(I) , BVY(I)
END DO
READ (1,*) NBT
DO I = 1 , NBT
READ (1,*) IBNDT (I) , BVT(I)
END DO
CLOSE (1)
C=======================================================================
C========> FILENAME XXXXXXX.BIN
OPEN ( 1, FILE=FNAME(2),STATUS='UNKNOWN',FORM='UNFORMATTED' )
READ (1) TIME, DT
READ (1) ( U(I) , I = 1 , NNODE )
READ (1) ( V(I) , I = 1 , NNODE )
READ (1) ( T(I) , I = 1 , NNODE )
CLOSE (1)
WRITE (*,*)' DT, DTMAX =', DT, DTMAX
WRITE (*,*)' TIME, TMAX =', TIME, TMAX
C=======================================================================
C========> FILENAME XXXXXXX.STM
OPEN ( 1, FILE=FNAME(3), STATUS='UNKNOWN' )
READ (1,*) NBS
READ (1,*) ( IBNDS(I), BVS(I) , I = 1 , NBS )
CLOSE (1)
C=======================================================================
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,RHS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* B(2,ND),RHS(MXN)
CHARACTER FILENAME*12
FILENAME = "STREAM.OUT"
CALL PLTLGO ( FILENAME )
CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,S,B,RHS)
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE TEMP ( MXE,MXN,ND,NE,NNODE,NODEX, XCOORD,YCOORD,
* T,S,B,TREF )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),
* B(2,ND),T(MXN)
CHARACTER FILENAME*12
FILENAME = "TEMP.OUT"
DO I = 1 , NNODE
T(I) = T(I) - TREF
END DO
CALL PLTLGO ( FILENAME )
CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,S,B,T)
DO I = 1 , NNODE
T(I) = T(I) + TREF
END DO
CALL PLTEXT
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
DO I = 1 , 4
X(I) = CRD(1,I)
Y(I) = CRD(2,I)
S(I) = SS(I)
END DO
DO LEVEL = 1 , NSTEP
SXY = START + (LEVEL-1) * DS
K = 1
DO I = 1 , 4
J = I + 1
IF ( I .EQ. 4 ) J = 1
IF ( (S(I)-SXY)*(S(J)-SXY) .LT. 0 ) THEN
T=(SXY-S(I))/(S(J)-S(I))
X0=X(J)*T+(1.D0-T)*X(I)
Y0=Y(J)*T+(1.D0-T)*Y(I)
IF ( K .GT. 0 ) THEN
CALL XMOVE ( X0, Y0 )
ELSE
CALL XDRAW ( X0, Y0 )
END IF
K = -K
END IF
END DO
END DO
RETURN
END
C
C
SUBROUTINE PLTUV (MXN,NNODE,XCOORD,YCOORD,U,V,NE,MXE,ND,NODEX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(MXN),YCOORD(MXN),U(MXN),V(MXN),NODEX(MXE,ND)
CHARACTER FILENAME*12
FILENAME = "VECTOR.OUT"
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 ( FILENAME )
CALL JCOLOR ( 11 )
CALL BOUND (MXE,MXN,ND,NE,NODEX,XCOORD, YCOORD)
CALL JCOLOR ( 11 )
CALL MINMAX ( MXN, NNODE, XCOORD, XMIN, XMAX )
CALL MINMAX ( MXN, NNODE, YCOORD, YMIN, YMAX )
FACT = RATIO * DMAX1 ( (YMAX-YMIN) , (XMAX-XMIN) ) / RMAX
DO I = 1, NNODE
CALL VECPLT ( XCOORD(I),YCOORD(I),U(I),V(I),FACT )
END DO
CALL PLTEXT
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 BOUND (MXE,MXN,ND,NE,NODEX,XCOORD, YCOORD)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(MXN),YCOORD(MXN), NODEX(MXE,ND)
LOGICAL LINE
C
DO IEL = 1 , NE
DO I = 1 , ND
LINE = .TRUE.
J = I + 1
IF ( I .EQ. ND ) J = 1
C
DO JEL = 1 , NE
IF ( IEL .NE. JEL ) THEN
DO II = 1 , ND
JJ = II + 1
IF ( II .EQ. ND ) JJ = 1
IF ( NODEX(IEL,I).EQ.NODEX(JEL,JJ) .AND.
* NODEX(IEL,J).EQ.NODEX(JEL,II) ) THEN
LINE = .FALSE.
EXIT
END IF
END DO
END IF
IF ( .NOT. LINE ) EXIT
END DO
IF ( LINE ) THEN
CALL XMOVE ( XCOORD(NODEX(IEL,I)) , YCOORD(NODEX(IEL,I)) )
CALL XDRAW ( XCOORD(NODEX(IEL,J)) , YCOORD(NODEX(IEL,J)) )
END IF
C
END DO
END DO
RETURN
END
C
C
SUBROUTINE PRESS (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* P,S,B )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),P(MXN),
* S(ND),B(2,ND)
CHARACTER FILENAME*12
FILENAME = "PRESSURE.OUT"
CALL PLTLGO ( FILENAME )
CALL CONTOUR ( MXE, MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,S,B,P)
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,YCOORD,
* S, B, P )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN),YCOORD(MXN),P(MXN),
* S(ND),B(2,ND)
C
CALL MINMAX ( MXN, NNODE, P, PPMIN, PPMAX )
IF ( PPMAX .EQ. PPMIN ) RETURN
NSTEP = 20
DS = ( PPMAX - PPMIN ) / NSTEP
I = PPMIN/DS+0.5D0
PPMIN = I*DS
C
CALL JCOLOR ( 9 )
IC = 0
CALL BOUND (MXE,MXN,ND,NE,NODEX,XCOORD, YCOORD)
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
CALL PLTSAI ( DS, NSTEP, PPMIN, B, S )
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
C
SUBROUTINE PLTEL4 ( ND,MXE,MXN,NODEX,NE,XCOORD, YCOORD )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(MXN), YCOORD(MXN)
CALL PLTLGO ("ELEMENT4.OUT")
DO I = 1 , NE
DO J = 1 , ND-1
CALL XMOVE ( XCOORD(NODEX(I,J )),YCOORD(NODEX(I,J )) )
CALL XDRAW ( XCOORD(NODEX(I,J+1)),YCOORD(NODEX(I,J+1)) )
END DO
CALL XMOVE ( XCOORD(NODEX(I,ND )),YCOORD(NODEX(I,ND )) )
CALL XDRAW ( XCOORD(NODEX(I, 1 )),YCOORD(NODEX(I, 1 )) )
END DO
CALL PLTEXT
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
C======================== GRAPHICS RUTINES =============================
SUBROUTINE PLTLGO ( FILENAME )
IMPLICIT REAL*8 ( A-H , O-Z )
CHARACTER FILENAME*12
OPEN ( 1, FILE=FILENAME, STATUS='UNKNOWN' )
RETURN
END
C
C
SUBROUTINE PLTEXT
CLOSE (1)
RETURN
END
C
C
SUBROUTINE JCOLOR ( I )
RETURN
END
C
C
SUBROUTINE XMOVE ( X , Y )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE(1,*) X, Y
RETURN
END
C
C
SUBROUTINE XDRAW ( X , Y )
IMPLICIT REAL*8 ( A-H , O-Z )
WRITE(1,*) X, Y
WRITE(1,*)
RETURN
END