PROGRAM NSG4FILE
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 25JAN. 2013
C=======================================================================
INCLUDE 'PARAM.DAT'
PARAMETER ( MXBA=MXB*2,ND2=ND*ND,INPT0=INTEPT-1 )
IMPLICIT REAL*8 ( A-H , O-Z )
C ARRAYS
CHARACTER*19 FNAME(NF), NAME
CHARACTER*15 PROJECT
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), IRPT(MXN)
DIMENSION IBNDS(MXB),BVS(MXB)
DIMENSION XCOORD(2,MXN),U(2,MXN),STM(MXN), P(MXN)
C=======================================================================
CALL GRULE ( INTEPT, SAI , W )
CALL GRULE ( INPT0 , SAI0, W0 )
CALL SHAPEF( ND, INTEPT, X, SAI , SF )
CALL DERIV ( ND, INTEPT, X, Y, SAI , BPP )
CALL DERIV ( ND, INPT0 , X, Y, SAI0, BP0 )
WRITE(*,*)' CFD GRAPHICS PROGRAM'
WRITE (*,*)' READING IN DATA FILES'
CALL INPUT ( MXB,MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,U,NF,
* FNAME,IBNDS,BVS,NBS,PROJECT,FLMDA )
WRITE (*,*)' PROJECT NAME =======>', PROJECT
C=======================================================================
CALL BANDWID ( MXE, ND, NE, NODEX, NBW )
IF ( NBW .GT. MXW ) STOP' ERROR #5'
CCC WRITE (*,*) 'COMPUTATION OF PRESSURE BY DOMAIN INTEGRATION'
CCC CALL COMPP ( MXE,MXN,INPT0,ND,BP0,NE, NNODE,FLMDA,BX,XCOORD,
CCC * U,NODEX,IRPT,P)
WRITE (*,*) 'COMPUTATION OF PRESSURE BY LINE INTEGRAL'
CALL COMPPGRN ( MXE,MXN,ND,NE,NNODE,FLMDA,XCOORD,U,NODEX,IRPT,P)
C=======================================================================
WRITE (*,*) 'COMPUTATION OF STREAMLINES'
CALL COMPSTM ( MXE,MXN,MXW,INTEPT,ND,BPP,W,SF,NE,NNODE,NBW,XCOORD,
* NODEX,SK,BX,U,AM,STM,MXB,NBS,BVS,IBNDS)
C=======================================================================
WRITE (*,*) 'CALL PLTEL4'
CALL PLTEL4 ( ND,MXE,MXN,NODEX,NE,XCOORD )
C=======================================================================
WRITE (*,*) 'CALL PLTUV'
CALL PLTUV ( MXN,NNODE,XCOORD,U,NE,MXE,ND,NODEX )
C=======================================================================
WRITE (*,*) 'CALL STREAM'
NSTEP = 25
CALL STREAM ( MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,BX,STM,NSTEP )
C=======================================================================
WRITE (*,*) 'CALL PRESS'
NSTEP = 30
CALL PRESS (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,BX,P,NSTEP )
C=======================================================================
STOP 'NORMAL 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 ( MXB,MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,U,NF,
* FNAME,IBNDS, BVS, NBS,PROJECT,FLMDA )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,MXN),U(2,MXN),
* IBNDS(MXB), BVS(MXB)
CHARACTER FNAME(NF)*19, PROJECT*15
LOGICAL YES
C------- RETURN VALUES: NEARLY ALL VARIBLES IN THIS SUBROUTINE
C========> FILENAME NSDATA.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)
WRITE (*,*) FNAME(I)
END DO
CLOSE (1)
WRITE(*,*)' PROJECT NAME:============>>>',PROJECT
C========> HEADER
OPEN ( 1, FILE=FNAME(1), STATUS='UNKNOWN' )
READ (1,*) VISCO, FLMDA
READ (1,*) TMAX
READ (1,*) MAXITE, ERMAX
C========> ELEMENT
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========> NODE
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(1,NODE) , XCOORD(2,NODE)
END DO
CLOSE (1)
C
C----------------------------------
C========> FILENAME XXXXXXX.BIN
OPEN ( 1, FILE=FNAME(2),STATUS='UNKNOWN',FORM='UNFORMATTED' )
READ (1) TIME, DT
READ (1) ( U(1,I) , I = 1 , NNODE )
READ (1) ( U(2,I) , I = 1 , NNODE )
CLOSE (1)
WRITE (*,*)' DT =', DT
WRITE (*,*)' TIME, TMAX =', TIME, TMAX
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)
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,NODEX,STIFF,B,U,STRM,RHS,MXB,NBS,BVS,IBNDS )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,MXN),
* SF(ND,INTEPT,INTEPT),STRM(MXN,MXW),STIFF(ND,ND),B(2,ND),U(2,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(1,NODE)
YAC12 = YAC12 + BPP(1,I,K,L) * XCOORD(2,NODE)
YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(1,NODE)
YAC22 = YAC22 + BPP(2,I,K,L) * XCOORD(2,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(1,NODE)
DVDX = DVDX + B(1,J) * U(2,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 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,U,NE,MXE,ND,NODEX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(2,MXN),U(2,MXN),NODEX(MXE,ND)
C-------------
CALL UMINMAX ( MXN, NNODE, XCOORD, XMIN, XMAX,YMIN,YMAX )
C-------------
CALL PLTLGO ( 'VECTOR00.OUT' )
PLTLMT = 0.001
SQLMT = PLTLMT **2
RATIO = 0.05
RMAX = 0.
DO I = 1 , NNODE
RMAX = DMAX1 ( RMAX, U(1,I)*U(1,I)+U(2,I)*U(2,I) )
END DO
IF ( RMAX .EQ. 0. ) RETURN
RMAX = DSQRT ( RMAX )
CALL UMINMAX ( MXN, NNODE, U, UMIN, UMAX,VMIN,VMAX )
CALL BOUND (MXE,MXN,ND,NE,NODEX,XCOORD)
FACT = RATIO * DMAX1 ( (YMAX-YMIN) , (XMAX-XMIN) ) / RMAX
DO I = NNODE , 1 , -1
IF ( (U(1,I)*U(1,I)+U(2,I)*U(2,I))/RMAX .GT. SQLMT )
* CALL VECPLT ( XCOORD(1,I),XCOORD(2,I),U(1,I),U(2,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 /
DX = U * FACT
DY = V * FACT
RNX = BETA * DY
RNY = - BETA * DX
WRITE (1,*) X0 , Y0
WRITE (1,*) X0+AL*DX , Y0+AL*DY
WRITE (1,*) X0+AL*DX+RNX , Y0+AL*DY+RNY
WRITE (1,*) X0+ DX , Y0+ DY
WRITE (1,*) X0+AL*DX-RNX , Y0+AL*DY-RNY
WRITE (1,*) X0+AL*DX , Y0+AL*DY
WRITE (1,*)
RETURN
END
C
C
SUBROUTINE STREAM (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,B,STM,NSTEP)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(2,MXN), B(2,ND),STM(MXN)
CHARACTER FILENAME*12
FILENAME = "STREAM00.OUT"
CALL PLTLGO ( FILENAME )
CALL CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,B,STM,NSTEP)
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE CONTOUR (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,B,P,NSTEP)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION S(ND),NODEX(MXE,ND),XCOORD(2,MXN),B(2,ND),P(MXN)
CALL MINMAX ( MXN, NNODE, P, PPMIN, PPMAX )
IF ( PPMAX .EQ. PPMIN ) RETURN
NSTEP = NSTEP/2*2 + 1
DS = ( PPMAX - PPMIN ) / NSTEP
I = PPMIN/DS
PPMIN= I*DS
C------------
CALL BOUND (MXE,MXN,ND,NE,NODEX,XCOORD)
DO IEL = 1 , NE
DO I = 1 , ND
B(1,I) = XCOORD(1,NODEX(IEL,I))
B(2,I) = XCOORD(2,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 PRESS (MXE,MXN,ND,NE,NNODE,NODEX,XCOORD,S,B,P,NSTEP )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,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,S,B,P,NSTEP)
CALL PLTEXT
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 UMINMAX ( MXN, NNODE, U, UMIN, UMAX,VMIN,VMAX )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION U(2,MXN)
UMIN = U(1,1)
UMAX = U(1,1)
VMIN = U(2,1)
VMAX = U(2,1)
DO I = 1 , NNODE
UMIN = DMIN1 ( UMIN , U(1,I) )
UMAX = DMAX1 ( UMAX , U(1,I) )
VMIN = DMIN1 ( VMIN , U(2,I) )
VMAX = DMAX1 ( VMAX , U(2,I) )
END DO
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,U,NODEX, IRPT, P )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,MXN), BX(2,ND),
* BPP(2,ND,INTEPT,INTEPT),U(2,MXN), IRPT(MXN),P(MXN)
C-------- COMPUTATION OF PRESSURE
C-------- RETURN VALUES P(MXN)
C-------- RESET
WRITE (*,*) FLMDA
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(1,NODE)
YAC12 = YAC12 + BPP(1,I,K,L) * XCOORD(2,NODE)
YAC21 = YAC21 + BPP(2,I,K,L) * XCOORD(1,NODE)
YAC22 = YAC22 + BPP(2,I,K,L) * XCOORD(2,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(1,NODE)+BX(2,J)*U(2,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
SUBROUTINE COMPPGRN ( MXE,MXN,ND,NE,NNODE,FLMDA,XCOORD,U,
* NODEX, IRPT, P )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,MXN),U(2,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---------- ELEMENT AREA COMPUTATION -----------
DO IEL = 1 , NE
ELEMENTA = 0.D0
DIVVINTE = 0.D0
DO ISEG = 1 , ND
J = ISEG + 1
IF ( ISEG .EQ. ND ) J = 1
I = NODEX(IEL,ISEG)
J = NODEX(IEL,J)
DX = XCOORD(1,J) - XCOORD(1,I)
DY = XCOORD(2,J) - XCOORD(2,I)
XM = 0.5D0*(XCOORD(1,I)+XCOORD(1,J))
YM = 0.5D0*(XCOORD(2,I)+XCOORD(2,J))
C
UM = 0.5D0*(U(1,I)+U(1,J))
VM = 0.5D0*(U(2,I)+U(2,J))
C
ELEMENTA = ELEMENTA + 0.5D0*(DY*XM-DX*YM)
DIVVINTE = DIVVINTE + (DY*UM-DX*VM)
END DO
C
DIVV = DIVVINTE / ELEMENTA
PRESS = - FLMDA*DIVV
C--------- DISTRIBUTION OF PRESSURE
DO J = 1 , ND
I = NODEX(IEL,J)
P(I) = P(I) + PRESS
IRPT(I) = IRPT(I) + 1
END DO
END DO
C
DO I = 1 , NNODE
P(I) = P(I) / IRPT(I)
END DO
RETURN
END
C
C
SUBROUTINE PLTEL4 ( ND,MXE,MXN,NODEX,NE,XCOORD )
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION NODEX(MXE,ND),XCOORD(2,MXN)
CALL PLTLGO ("ELEMENT4.OUT")
DO I = 1 , NE
DO J = 1 , ND-1
CALL XMOVE ( XCOORD(1,NODEX(I,J )),XCOORD(2,NODEX(I,J )) )
CALL XDRAW ( XCOORD(1,NODEX(I,J+1)),XCOORD(2,NODEX(I,J+1)) )
END DO
CALL XMOVE ( XCOORD(1,NODEX(I,ND )),XCOORD(2,NODEX(I,ND )) )
CALL XDRAW ( XCOORD(1,NODEX(I, 1 )),XCOORD(2,NODEX(I, 1 )) )
END DO
CALL PLTEXT
RETURN
END
C
C
SUBROUTINE BOUND (MXE,MXN,ND,NE,NODEX,XCOORD)
IMPLICIT REAL*8 ( A-H , O-Z )
DIMENSION XCOORD(2,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(1,NODEX(IEL,I)) , XCOORD(2,NODEX(IEL,I)) )
CALL XDRAW ( XCOORD(1,NODEX(IEL,J)) , XCOORD(2,NODEX(IEL,J)) )
END IF
C
END DO
END DO
RETURN
END
C
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
IMPLICIT REAL*8 ( A-H , O-Z )
CLOSE (1)
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